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Foreword 


There are situations in clinical reality when it would be beneficial to be able to use a 
structural and functional prosthesis to compensate for a congenital or acquired 
defect that can not be replaced by biologic material. 

Mechanical stability of the connection between material and biology is a 
prerequisite for successful rehabilitation with the expectation of life long function 
without major problems. 

Based on Professor Skalak’s theoretical deductions of elastic deformation at/of 
the interface between a screw shaped element of pure titanium at the sub cellular 
level the procedure of osseointegration was experimentally and clinically developed 
and evaluated in the early nineteen-sixties. 

More than four decades of clinical testing has ascertained the predictability of 
this treatment modality, provided the basic requirements on precision in 
components and procedures were respected and patients continuously followed. 

The functional combination of a piece of metal with the human body and its 
immuno-biologic control mechanism is in itself an apparent impossibility. Within 
the carefully identified limits of biologic acceptability it can however be applied 
both in the cranio-maxillofacial skeletal as well as in long bones. 

This book provides an important contribution to clinical safety when bone 
anchored prostheses are used because it explains the mechanism and safety margins 
of transfer of load at the interface with emphasis on the actual clinical anatomical 
situation. This makes it particularly useful for the creative clinician and unique in its 
field. It should also initiates some critical thinking among hard ware producers who 
might sometimes underestimate the short distance between function and failure when 
changes in clinical devices or procedures are too abruptly introduced. 

An additional value of this book is that it emphasises the necessity of respect 
for what happens at the functional interaction at the interface between molecular 
biology and technology based on critical scientific exploration and deduction. 


P-I Branemark 


Preface 


This book provides the theoretical foundation of Finite Element Analysis(FEA) in 
implant dentistry and practical modelling skills that enable the new users (implant 
dentists and designers) to successfully carry out FEA in actual clinical situations. 
The text is divided into five parts: introduction of finite element analysis and 
implant dentistry, applications, theory with modelling and use of commercial 
software for the finite element analysis. The first part introduces the background of 
FEA to the dentist in a simple style. The second part introduces the basic 
knowledge of implant dentistry that will help the engineering designers have some 
backgrounds in this area. The third part is a collection of dental implant applications 
and critical issues of using FEA in dental implants, including bone-imp lant interface, 
implant-prosthesis connection, and multiple implant prostheses. The fourth part 
concerns dental implant modelling, such as the assumptions of detailed geometry of 
bone and implant, material properties, boundary conditions, and the interface 
between bone and implant. Finally, in fifth part, two popular commercial finite 
element software ANSYS and ABAQUS are introduced for a Branemark same-day 
dental implant and a GJP biomechanical optimum dental implant, respectively. 


Jianping Geng Hangzhou 
Weiqi Yan Hangzhou 
Wei Xu Surrey 
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4.1 Introduction 


The use of numerical methods such as FEA has been adopted in solving complicated 
geometric problems, for which it is very difficult to achieve an analytical solution. 
FEA is a technique for obtaining a solution to a complex mechanics problem by 
dividing the problem domain into a collection of much smaller and simpler domains 
(elements) where field variables can be interpolated using shape functions. An 
overall approximated solution to the original problem is determined based on 
variational principles. In other words, FEA is a method whereby, instead of seeking 
a solution function for the entire domain, it formulates solution functions for each 
finite element and combines them properly to obtain a solution to the whole body. 
A mesh is needed in FEA to divide the whole domain into small elements. The 
process of creating the mesh, elements, their respective nodes, and defining 
boundary conditions is termed “discretization” of the problem domain. Since the 
components in a dental implant-bone system is an extremely complex geometry, 
FEA has been viewed as the most suitable tool to mathematically, model it by 
numerous scholars. 
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FEA was initially developed in the early 1960s to solve structural problems in 
the aerospace industry but has since been extended to solve problems in heat 
transfer, fluid flow, mass transport, and electromagnetic realm. In 1977, Weinstein! 
was the first to use FEA in implant dentistry. Subsequently, FEA was rapidly 
applied in many aspects of implant dentistry. Atmaram and Mohammed** analysed 
the stress distribution in a single tooth implant, to understand the effect of elastic 
parameters and geometry of the implant, implant length variation, and pseudo- 
periodontal ligament incorporation. Borchers and Reichart’ performed a three 
dimensional FEA of an implant at different stages of bone interface development. 
Cook, et al.° applied it in porous rooted dental implants. Meroueh, et al.’ used it for 
an osseointegrated cylindrical implant. Williams, et al.° carried out it on cantilevered 
prostheses on dental implants. Akpinar, et al.’ simulated the combination of a nature 
tooth with an implant using FEA. 


4.2 Considerations of Dental Implant FEA 


In the past three decades, FEA has become an increasingly useful tool for the 
prediction of stress effect on the implant and its surrounding bone. Vertical and 
transverse loads from mastication induce axial forces and bending moments and 
result in stress gradients in the implant as well as in the bone. A key to the success 
or failure of a dental implant is the manner in which stresses are transferred to the 
surrounding bone. Load transfer from the implant to its surrounding bone depends 
on the type of loading, the bone-implant interface, the length and diameter of the 
implants, the shape and characteristics of the implant surface, the prosthesis type, 
and the quantity and quality of the surrounding bone. FEA allows researchers to 
predict stress distribution in the contact area of the implant with cortical bone and 
around the apex of the implant in trabecular bone. 

Although the precise mechanisms are not fully understood, it is clear that there 
is an adaptive remodelling response of the surrounding bone to this kind of stress. 
Implant features causing excessive high or low stresses can possibly contribute to 
pathologic bone resorption or bone atrophy. The principal difficulty in simulating 
the mechanical behaviour of dental implants is the modelling of human bone tissues 
and its response to applied mechanical forces. The complexity of the mechanical 
characterization of bone and its interaction with implant systems have forced 
researchers to make major simplifications and assumptions to make the modelling 
and solving process possible. Some assumptions influence the accuracy of the FEA 
results significantly. They are: (1) detailed geometry of the bone and implant to be 
modelled’®, (2) material properties’®, (3) boundary conditions’’, and (4) the interface 
between the bone and implant”’. 
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4.3 Fundamentals of Dental Implant Biomechanics 


4.3.1 Assumptions of Detailed Geometry of Bone and Implant 


The first step in FEA modelling is to represent the geometry of interest in the 
computer. In some two- or three-dimensional FEA studies the bone was modeled as 
a simplified rectangular configuration with the implant'’” (Fig. 4.1). Some three- 
dimensional FEA models treated the mandible as an arch with rectangular section’*”’. 
Recently, with the development of digital imaging techniques, more efficient 
methods are available for the development of anatomically accurate models. These 
include the application of specialized softwares for the direct transformation of 2D 
or 3D information in image data from CT or MRI, into FEA meshes (Fig.4.2 to Fig. 
4.4). The automated inclusion of some material properties from measured bone 
density values is also possible’*’’. This will allow more precise modelling of the 
geometry of the bone-implant system. In the foreseeable future, the creation of FEA 
models for individual patients based on advanced digital techniques will become 
possible and even commonplace. 


Abutment 


——. Abutment screw 


Fixture body 


INNA AN 


Fig. 4.1 3D Information of a Simplified Rectangular Configuration with the Implant 
Components (By H.M. Huang and S.Y. Lee) 
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Fig. 4. 2 2D Information in Image Data Fig. 4.3 3D Information in 
from Mandibular CT and Its FEA Stress Image Data from Posterior Maxillary 


Distribution (By J.P. Geng, et al.) CT and Its FEA Meshes 
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(a) (b) 


Fig. 4.4 3D Information in Image Data from Mandibular CT and Its FEA Meshes 


4.3, 2 Material Properties 


Material properties greatly influence the stress and strain distribution in a structure. 
These properties can be modeled in FEA as isotropic, transversely isotropic, 
orthotropic, and anisotropic. In an isotropic material, the properties are the same in 
all directions, and therefore there are only two independent material constants. An 
anisotropic material has its different properties when measured in different 
directions. There are many material constants depending on the degree of anisotropy 
(transversely isotropic, orthotropic, etc). 

In most reported studies, materials are assumed homogeneous, linear and have 
elastic material behaviour characterized by two material constants of Young's 
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modulus and Poisson’s Ratio. Early FEA studies ignored the trabecular bone 
network simply because it’s pattern was not able to be determined. Therefore, it 
was assumed that trabecular bone has a solid pattern inside the inner cortical bone 
shell. Both bone types were simplistically modeled as linear, homogeneous, and 
isotropic materials. A range of different material parameters have been recommended 
for use in previous FEA studies (Table 4.1)". 


Table 4. 1 Material Parameters Used in FEA Studies of Dental Implants 


Ap Mois PoSOn'S Author, Year ———_—Reference No. 
Enamel 4.14 X10! 0.3 Davy, 1981 18 
4.689 X 10° 0.30 Wright, 1979 19 
8.25 X 10° 0.33 Farah, 1975 20 
8.4 X10 0.33 Farah, 1989 21 
Dentin 1.86 X10 0.31 Reinhardt, 1984 22 
1.8 X 10° 0.31 MacGregor, 1980 23 
Parodontal 17] 0.45 Atmaram, 198] 24 
membrance 69.8 0.45 Reinhardt, 1984 22 
6.9 0.45 Farah, 1989 21 
Cortical bone 2727 0.30 Rice, 1988 25 
1.0 X 10° 0.30 Farah, 1989 21 
1.34 X10 0.30 Cook, 1982 26 
1.5 10° 0.30 Cowin, 1989 27 
Trabecular bone 150 0.30 Cowin, 1989 27 
250 0.30 MacGregor, 1980 23 
790 0.30 Knoell, 1977 28 
1.37 X10 0.31 Borchers, 1983 5 
Mucosa 10 0.40 Maeda, 1989 29 
Pure Ti 117 X10° 0.30 Ronald, 1995 30 
Tr6Al4V 110 X10 0.33 Colling, 1984 31 
Type 3 gold alloy 100 X 10° 0.30 Ronald, 1995 30 
80 X10 0.33 Lewinstein, 1995 32 
Ag-Pd alloy 95 X10 0.33 Craig, 1989 33 
Co-Cr alloy 218 X10 0.33 Craig, 1989 33 
Porcelain 68.9 X10 0.28 Lewinstein, 1995 32 
Resin 29 AAG 0.35 Craig, 1989 33 
Resin composite 7 X10 0.2 Craig, 1989 33 


In fact, several studies**’ have pointed out that cortical bone is neither 
homogeneous nor isotropic (Table 4. 2). This non-homogenous, anisotropic, 
composite structure of bone also possesses different values for ultimate strain and 
modulus of elasticity when it 1s tested in compression compared to in tension. Test 
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conditions will affect the material properties measured too. Rieger, et al.” reported 
that a range of stresses (1.4 MPa to 5.0 MPa) appeared to be necessary for healthy 
maintenance of bone. Stresses outside this range have been reported to cause bone 
resorption. 


Table 4, 2 Anisotropic Properties of Cortical Bone 


Cortical Shell 
Elastic (M Pa) 


Diap hy seal Metaphyseal 
Longitudinal 17,500 9,650 
Transverse 11,500 5,470 


4.3.3 Boundary Conditions 


Most FEA studies modelling the mandible set boundary conditions as a fixed 
boundary. Zhou” developed a more realistic three-dimensional mandibular FEA 
model from transversely scanned CT image data. The functions of the muscles of 
mastication and the ligamenteous and functional movement of the TMJs were 
simulated by means of cable elements and compressive gap elements respectively. It 
was concluded that cable and gap elements could be used to set boundary conditions 
in their mandibular FEA model, improving the model mimicry and accuracy. 
Expanding the domain of the model can reduce the influence of inaccurate modelling 
of the boundary conditions. This however, will be at the expense of computing and 
modelling time (Fig. 4.5). Teixera, et al.” concluded that in a three-dimensional 
mandibular model, modelling the mandible at distances greater that 4. 2 mm mesial or 
distal from the implant did not result in any significant further yield in FEA 
accuracy. Use of infinite elements can be a good way to model boundary conditions 


(Fig.4.6). 


Fig. 4. 5 Three-dimensional FEA Model of the Human Jaw and the Functional 
Direction of the Muscles of Mastication 


4.4 Interface between Bone and Implant 


FEA models usually assume a state of optimal osseointegration, meaning that 
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Fig. 4.6 Illustration of the Base Model of the Mandible (left) and the Longest Model 
(right) (Reproduced from J Oral Rehabil 1998;25:300 with permission) 


cortical and trabecular bone is perfectly bonded to the implant. This does not occur 
exactly in clinical situation. Therefore, the imperfect contact and its effect on load 
transfer from implant to supporting bone need to be modelled more carefully. 
Current FEA programs provide several types of contact algorithms to technically 
conduct simulation of clinical contacts. The friction between contact surfaces can 
also be modelled with contact algorithms. The friction coefficients, however, have to 
be determined via experimentation. 

Bone is a porous material with complex microstructure. The higher load bearing 
capacity of dense cortical bone compared to the more porous trabecular bone is 
generally recognised. Upon implant insertion, cortical and/or trabecular bone, starting 
at the periosteal and endosteal surfaces, gradually forms a partial to complete 
encasement around the implant. However, the degree of encasement is dependent on 
the stresses generated and the location of the implant in the jaw’. The anterior 
mandible is associated with 100% cortical osseointegration and this decreases 
toward the posterior mandible. The least cortical osseointegration (< 25% ) is seen 
in the posterior maxilla. The degree of osseointegration appears to be dependant on 
bone quality and stresses developed during healing and function. To study the 
influence of osseointegration in greater detail at the bone trabeculae contact to 
implant level, Sato, et al.” set up four types of stepwise assignment algorithms of 
elastic modulus according to the bone volume in the cubic cell (Fig.4.7 to Fig.4.8). 
They showed that a 300 um element size was valid for modelling the bone-imp lant 
interface. 
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Fig. 4.7 Construction Procedure of Bone Trabecular Model (Reproduced from J Oral 
Rehabil 1998;26:641 with permission) 
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Fig. 4. 8 Four Types of Stepwise Assignment Algorithms of Young’s Modulus 
According to the Bone Volume in the Cubic Cell (E: elastic modulus, GPa) (Reproduced from 
J Oral Rehabil 1999;26:641 with permission) 


4.5 Reliability of Dental Implant FEA 


Stress distribution depends on assumptions made in modelling geometry, material 
properties, boundary conditions, and bone-implant interface. To obtain more 
accurate stress predictions, advanced digital imaging techniques can be applied in 
modelling more realistical bone geometry; the anisotropic and nonhomogenous nature 
of materials need to be considered; and boundary conditions have to be carefully 
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treated using computational modelling techniques. In addition, modelling of the bone- 
implant interface should incorporate the actual osseointegration contact area in 
cortical bone as well as the detailed trabecular bone contact pattern, using contact 
algorithms in FEA. 


4.6 Conclusions 


FEA has been used extensively in the prediction of biomechanical performance of 
dental implant systems. Assumptions made in the use of FEA in Implant Dentistry 
have to be taken into account when interpreting the results. 

FEA is an effective computational tool that has been used to analyse dental 
implant biomechanics. Many optimisations of design feature have been predicted 
and will be applied to potential new implant systems in the future. 
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3.1 Introduction 


Although the precise mechanisms are not fully understood, it is clear that there is an 
adaptive remodelling response of surrounding bone to stresses. Implant features 
causing excessive high or low stresses can possibly contribute to pathologic bone 
resorption or bone atrophy. This chapter reviews the current applications of FEA 
in Implant Dentistry. Findings from FEA studies will then be discussed in relation 
to the bone-implant interface; the implant-prosthesis connection; and multiple 
implant prostheses. 


3.2 Bone-implant Interface 


3.2.1 Introduction 


Analyzing force transfer at the bone-implant interface is an essential step in the 
overall analysis of loading, which determines the success or failure of an implant. It 
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has long been recognized that both implant and bone should be stressed within a 
certain range for physiological homeostasis. Overload can cause bone resorption or 
fatigue failure of the implant whilst underload may lead to disuse atrophy and to 
subsequent bone loss as well'’. Using load cells in rabbit calvaria, Hassler, et al.’ 
showed that the target compressive stress level for maximum bone growth occurs at 
1.8 MPa leveling off to a control level at 2.8 MPa. Skalak* states that close 
apposition of bone to the titantum implant surface means that under loading, the 
interface moves as a unit without any relative motion and this is essential for the 
transmission of stress from the implant to the bone at all parts of the interface. 

In centric loading, several FEA studies’'’ of osseointegrated implants 
demonstrated that when maximum stress concentration is located in cortical bone, it 
is in the contact area with the implant; and when maximum stress concentration is in 
trabecular bone, it occurs around the apex of the implant. In cortical bone, stress 
dissipation is restricted to the immediate surroundings of the implant, whereas in 
trabecular bone a fairly broader distant stress distribution occurs”. 


3.2.2 Stress Transmission and Biomechanical Implant Design Problem 


FEA can simulate the interaction phenomena between implants and the surrounding 
tissues. Analysis of the functional adaptation process is facilitated by accessing 
various loadings and implant and surrounding tissue variables. Load transfer at the 
bone-implant interface depends on (1) the type of loading, (2) material properties of 
the implant and prosthesis, (3) implant geometry-length, diameter as well as shape, 
(4) implant surface structure, (5) the nature of the bone-implant interface, and (6) 
quality and quantity of the surrounding bone. Most efforts have been directed at 
optimizing implant geometry to maintain the beneficial stress level in a variety of 
loading scenarios. 


3,2,2. 1 Loading 


When applying FEA to dental implants, it is important to consider not only axial 
forces and horizontal forces (moment-causing loads), but also a combined load 
(oblique bite force), since these are more realistic bite directions and for a given force 
will cause the highest localized stress in cortical bone''. Barbier, et al.” investigated 
the influence of axial and non-axial occlusal loads on the bone remodelling around 
IMZ implants in a dog mandible simulated with FEA. Strong correlation between 
the calculated stress distributions in the surrounding bone tissue and the remodelling 
phenomena in the comparative animal model was observed. They concluded that the 
highest bone remodelling events coincide with the regions of highest equivalent 
stress and that the major remodelling differences between axial and non-axial loading 
are largely determined by the horizontal stress component of the engendered 
stresses. The importance of avoiding or minimizing horizontal loads should be 
emp hasized. 

Zhang and Chen’ compared dynamic with static loading, in three-dimensional 
FEA models with a range of different elastic moduli for the implant. Their results 
showed that compared to the static load models, the dynamic load model resulted in 
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higher maximum stress at the bone-implant interface as well as a greater effect on 
stress levels when elastic modulus was varied. 

In summary, both static and dynamic loading of implants have been modelled 
with FEA. In static load studies, it is necessary to include oblique bite forces to 
achieve more realistic modelling. Most studies concluded that excessive horizontal 
force should be avoided. The effects of dynamic loading requires further 
investigation. 
3.2.2.2 Material Properties 


Prosthesis material properties 
High rigidity prostheses are recommended because the use of low elastic moduli 
alloys for the superstructure predicts larger stresses at the bone-implant interface 
on the loading side than the use of a rigid alloy with the same geometry’. Stegariou, 
et al.” used three-dimensional FEA to assess stress distribution in bone, implant, 
and abutment when gold alloy, porcelain, or resin (acrylic or composite) was used 
for a 3-unit prosthesis. In almost all cases, stress at the bone-implant interface with 
the resin prostheses was similar to, or higher than that in the models with the other 
two prosthetic materials. But in his classical mechanical analysis, Skalak’* stated 
that the presence of a resilient element in an implant prosthesis superstructure 
would reduce the high load rates that occur when biting unexpectedly on a hard 
object. For this reason, he suggested the use of acrylic resin teeth. Nevertheless, 
several other studies*’’ could not demonstrate any significant differences in the 
force absorption quotient of gold, porcelain or resin prostheses. 
Implant material properties 
The elastic moduli of different implant materials will influence the implant-bone 
interface. Implant materials with too low moduli should be avoided and Malaith, et 
al.'* suggested implant materials with an elastic modulus of at least 110,000 N/mm’. 
Rieger, et al.” indicated that serrated geometry led to high-stress concentrations at 
the tips of the bony ingrowth and near the neck of the implant. Low moduli of 
elasticity strengthened these concentrations. The non-tapered screw-type geometry 
showed high-stress concentrations at the base of the implant when high moduli were 
modeled and at the neck of the implant when low moduli were modeled. The authors 
concluded that a tapered endosseous implant with a high elastic modulus would be 
the most suitable for dental implantology. However, the design must not cause high- 
stress concentrations at the implant neck that commonly leads to bone resorption. 
Stoiber” reported that in the construction of an appropriate screw implant, special 
attention must be paid to high rigidity of the implant, rather than to thread design. 
In summary, although the effect of prosthesis material properties is still under 
debate, it has been well established that implant material properties can greatly 
affect the location of stress concentrations at the implant-bone interface. 


3.2.2.3 Implant Geometry 


Implant diameter and length 
Large implant diameters provide for more favourable stress distributions”. FEA 
has been used to show that stresses in cortical bone decrease in inverse proportion 
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to its increase in implant diameter with both vertical and lateral loads”. However, 
Holmgren, et al.’ showed that using the widest diameter implant is not necessarily 
the best choice when considering stress distribution to surrounding bone, but within 
certain morphological limits, an optimum dental implant size exists for decreasing 
the stress magnitudes at the bone-implant interface. 


(b) 


Fig, 3.1 GJP Implant Shape and Its Thread Design(By J.P. Geng) 


In general, the use of short implants has not been recommended because it is 
believed that occlusal forces must be dissipated over a large implant area to preserve 
the bone. Lum” has shown that occlusal forces are distributed primarily to the 
crestal bone rather than evenly throughout the entire surface area of the implant 
interface. Since masticatory forces are light and fleeting, these forces are normally 
well-tolerated by the bone. It is the bruxing forces that must be adequately 
attenuated, and this may be done by increasing the diameter and number of implants. 
A recent clinical study concluded that short implants are possible when the pert 
implant tissues were in good condition”. 

In summary, the optimum length and diameter necessary for long term 
implantation success depends on the bone support condition. If the bone is in 
normal condition, length and diameter appear not to be significant factors for 
implant success. However, if the bone condition is poor, large diameter implants are 
recommended and short implants should be avoided. 

Implant shape 

Holmgren, et al."' report that a stepped cylindrical design for press-fit situations is 
most desirable from the standpoint of stress distribution to surrounding bone. Using 
FEA to analyse a parasaggital model digitized from a computed tomography (CT) 
generated patient data set, these authors simulated various single-tooth, two- 
dimensional osseointegrated dental implant models. The results suggested that stress 
is more evenly dissipated throughout the stepped cylindrical implant compared to 
the straight implant type. Rieger, et al.’ concluded that a tapered endosseous 
implant with a high elastic modulus would be the most suitable one after analysing 
stress concentration patterns using FEA. Also using FEA, Mailath, et al.’* compared 
cylindrical to conical implant shapes when exposed to physiologic stresses and 
examined the occurrence of stress concentrations at the site of implant entry into 
bone. They reported that cylindrical implants were preferable to conical shapes. 
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Siegele and Soltesz“ compared cylindrical, conical, stepped, screw and hollow 
cylindrical implant shapes by means of FEA. Both a fixed bond (simulating 
complete load transfer with bioactive materials) and a pure contact (only 
compression transfer with biomert materials) without friction between implant and 
bone were considered as interface conditions. The results demonstrate that different 
implant shapes lead to significant variations in stress distributions in the bone. The 
authors stated that implant surfaces with very small radii of curvature (conical) or 
geometric discontinuities (stepped) induced distinctly higher stresses than smoother 
shapes (cylindrical, screw-shaped). Moreover, a fixed bond between implant and 
bone in the medullary region (as may be obtained with a bioactive coating) is 
advantageous to the stress delivered to bone, since it produces a more uniform stress 
distribution than a pure contact does. Clift, et al.” reported that the modification of 
the standard implant design to include a flexible central post resulted in a decrease in 
the maximum von Mises stresses and equivalent strains in cancellous bone. It was 
postulated that this would reduce the likelihood of bone fatigue failure and 
subsequent bone resorption. 

Optimum implant shape is related to the bone condition and implant material 
properties. Implant designs have adopted various shapes and FEA seems to indicate 
that for commercially pure titanium (CPTi) implants, smoother profiles engender 
lower stress concentrations. The optimal thread design to achieve the best load 
transfer characteristics is the subject of current investigations. 


3.2.2.4 Implant Surface Structure 


(500x] 500x [500x] 


Fig. 3.2 Etched Titanium, Hydroxylapatite-coated (HA-coated) Titanium, and Titanium 
Plasma Spray (TPS) Surface Treatments (Steri-Oss Implant System) 


Bioactive materials are used as coating on titanium implants because they have the 
potential to encourage bone growth up to the surface of the implant”. It is claimed 
that these coatings can produce a fully integrated interface with direct bonding 
between bone and the implant material, leading to a more even transfer of load to the 
bone along the implant and thus a reduction in stress concentrations”. 

Polyactive is a system of poly (ethylene oxide) poly (butylene terephthalate) 
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segmented co-polymers with bone-bonding capacity. Meijer, et al.”” investigated the 
influence of a three-layered flexible coating of Polyactive on bone stress distribution 
using three-dimensional FEA in a mandibular model. In the case of sagittal and 
transversal loading, the use of a Polyactive coating reduced both the minimum 
principal stress in the bone and the compressive radial stress at the bone-implant 
interface. However, it raised the maximum principal and the tensile radial stress. In 
the case of vertical loading, the application of a flexible coating reduced the 
compressive radial stress at the bone-implant interface around the neck of the 
implant by a factor of 6.6 and the tensile radial stress by a factor of 3.6. Variations 
in composition and thickness of the coating did not significantly affect the results. 


3.2.2.5 Nature of Bone-implant Interface 


There are two types of contact at the bone-implant interface: bone/implant contact 
and fibrous tissue/mplant contact. The clinical concept of fibrous encapsulation of 
an implant is considered to be a failure and this condition is no longer modelled in 
FEA studies. 


3.2.2.6 Surrounding Bone Quality and Quantity 


The long-term clinical performance of a dental implant is dependent upon the 
preservation of good quality bone surrounding the implant and a sound interface 
between the bone and the biomaterial. Good quality bone itself relies on an 
appropriate level of bone remodelling necessary to maintain the bone density and 
the avoidance of bone microfracture and failure. Both processes are governed by the 
stress and strain distribution in the bone. 

Crestal bone loss 

The crestal bone region is of particular interest due to observation of progressive 
bone resorption (saucerization). Crestal loss is observed around various designs of 
dental implants. A possible cause of this loss is related to the low stresses acting on 
pertimplant bone. An equivalent stress of 1.6 MPa is determined to be sufficient to 
avoid crestal bone loss from disuse atrophy in the canine mandibular premolar 
region, based on both histological examination and FEA results”*”” 

Wiskott and Belser” studied the relationship between the stresses applied and 
bone homeostasis of different implant neck designs. They observed that polished 
necks of dental implants did not osseointegrate as do textured surfaces. Lack of 
osseointegration was postulated to be due to increased pressure on the osseous bed 
during implant placement, establishment of a physiological “biologic width”, stress 
shielding, and lack of adequate biomechanical coupling between the load-bearing 
implant surface and the surrounding bone. Any viable osseous structure (including 
the tissue that surrounds the polished implant neck) is subjected to periodic phases 
of resorption and formation. Hansson’! compared implants with smooth necks to 
implants with retention elements all the way up to the crest. His FEA study found 
that retention elements at the implant neck resulted in major decrease in peak 
interfacial shear stresses. He suggests that these retention elements at the implant 
neck will counteract marginal bone resorption in accordance with Wolff's law. 

For the Screw-Vent implant, Clelland, et al.” showed that under axial loading, 
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mesial and distal stresses were much lower than those buccal and lingual to the 
implant. Maximum stress in the bone was lingual to the superior portion of the 
collar. Previous longitudinal radiographic studies of a similar implant had revealed 
bone loss mesial and distal to the implant. The authors conclude that the clinical 
significance of the stress transfer to the bone buccal and lingual to the implant has 
yet to be determined. 

Minimum required load for avoidance of crestal bone loss appears to have been 
defined *”, but the upper limit of the physiological stress range has not yet been 
investigated. 

Cortical bone 

The quality and quantity of the surrounding bone influences the load transfer from 
implant to bone’*’. In almost all FEA studies of titanium implants, stress 
concentrations occur around the implant neck. Under oblique loads with high 
occlusal stress magnitudes, the elastic limit of bone surrounding implants may be 
surpassed and lead to microfractures in the cortical bone. Clift, et al.’ emphasized 
the importance of having good quality dense bone around the implant neck which 
can withstand stresses in the range of 9-18 MPa prior to loading. Failure to achieve 
this after implantation and subsequent healing may result in local fatigue failure and 
resorption at the neck upon resumption of physiological loading*. Holmes and 
Loftus” examined the influence of bone quality on the transmission of occlusal 
forces for endosseous dental implants using FEA. Placement of implants in bone 
with greater thickness of the cortical shell and greater density of the core will result 
in less micro-movement and reduced stress concentration, thereby imcreasing the 
likelihood of fixture stabilization and tissue integration. 

Papavasiliou, et al.° showed with a three-dimensional FEA model that the 
absence of cortical bone increased interfacial stresses at the locations studied. Clift, 
et al.“ reported that a reduction in the elastic modulus of the bone around the neck 
of the implant by a factor of 16 only produced a two-fold reduction in the peak 
stress. 

Trabecular (cancellous) bone 

Using the degree of direct bone-implant interface as an indicator of endosseous 
implant success appears to have been over interpreted because 100% bone 
apposition is almost never obtained at the surface of the endosseous dental implant. 
Investigating the three-dimensional bone interface to hydroxy apatite-coated titanium 
alloy implants, Wadamoto, et al.’ generated computer graphics by the integration of 
data for serial ground surfaces obtained at 75 um intervals of the tissue block 
involved with the implant. They found that the bone contact ratio of the whole 
surface of each of the three implants was 80.8% , 68.1% , and 68.8% , and that for 
each direction and portion varied with the conditions of implant placement. The 
bone volume ratios around the implant at the 0 to 300 4m zone were also calculated, 
and total ratios ranged from 58% to 81%. These results may provide useful 
quantitative information about the bone structure around implants and contribute to 
the development of more realistic FEA models based on the biologic bone structure 
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around implants. 

Clelland, et al.* investigated a Steri-Oss implant in various bone models with 
different cancellous and cortical bone conditions using two-dimensional FEA. For 
the all-cancellous bone model, low stresses and high strains surrounded the implant 
apex. For the models with a layer of cortical bone added, higher crestal stresses and 
lower apical strains were observed. The thick layer(3mm) of isotropic cortical bone 
produced stresses at least 50% less than the thin layer(1. 5mm). The assumption of 
transverse isotropy (orthotropy) for cortical bone layer increased stresses and 
strains by approximately 25% compared with isotropic bone. They conclude that 
crestal cortical layer thickness and bone isotropy have a substantial impact on 
resultant stresses and strains. 


3.2.3 Summary 


Load transmission and resultant stress distribution is significant in determining the 
success or failure of an implant. Factors that influence the load transfer at the bone- 
implant interface include the type of loading, material properties of implant and 
prosthesis, length and diameter of implants, implant shape, structure of implant 
surface, nature of bone-implant interface, and the quality and quantity of 
surrounding bone. Of these biomechanical factors, implant length, diameter, and 
shape are easily changed. The quality and quantity of cortical and cancellous bone 
needs to be assessed clinically and will influence implant selection. 


3.3 Implant Prosthesis Connection 


3.3.1 Introduction 


Clinical studies have reported a significant incidence of component failure. These 
include gold screw failures, abutment screw failures, gold cylinder fractures, 
framework fractures, and implant fractures. The etiology of these failures is complex 
and involves cyclic fatigue, oral fluids, and varied chewing patterns and loads. 
Biomechanically, the following component interfaces can be found in the 
Branemark implant: (1) ftxture-abutment interface; (2) abutment screw-abutment 
interface; (3) gold cylinder-abutment interface; (4) gold retaining screw-gold cylinder 
interface; (5) gold retaining screw-abutment screw interface. Long term screw joint 
integrity at the fixture-abutment screw joint and abutment-gold cylinder screw joint 
is essential for prosthetic success. An increasing number of FEA studies are 
focusing on biomechanical problems involving the screw joint and on screw 


loosening phenomena”. 


3.3.2 Screw Loosening Problem 


3.3.2.1 Introduction 


The screw loosening problem frequently affects dental implants and imp lant- 
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supported prostheses. When a screw is fastened to fix a prosthesis, a tensile force 
(preload) is built up in the shank of the screw. This preload acts on the screw shank 
from the head of the screw to the threads. It should be as high as possible because it 
creates a clamping force between the abutment and the implant. The screw elongates 
when subjected to tensile forces during tightening. The more elongation there is, the 
better the stability of the screw in place. Thus screw design is of significance and 
should allow a maximum torque to be introduced on the shank of the screw”. 

Several authors’” have drawn attention to the fact that repeated loading and 
unloading cycles result in alternating contact and separation of components. Clinical 
findings of screw loosening and failure probably result from these separation events 
and from elevated strains in the screw. Another mechanism of screw loosening is 
based on the fact that no surface is completely smooth. Because of the 
microroughness of component surfaces, when the screw interface is subjected to 
external loads, micromovements occur between the surfaces. Wear of the contact 
areas might result from these motions, thereby bringing the two surfaces closer to 
each other and cause a decrease in preload in the set of screws. 

With prosthesis superstructure distortion, an external preload can be 
superimposed on the screw joints of the implant prosthesis. This distortion (or lack 
of passive fit) can impart additional axial forces and bending moments on the screw 
joints and increase the like'ihood of prosthetic component failure’. 


3.3.2.2 Application of Preload 


The load-transfer mechanism between prosthetic components arises from torque 
application to the abutment screw and gold screw. Sakaguchi, et al.” developed a 
two-dimensional FEA model for non-linear contact analysis of Branemark implant 
prosthetic components. They found that when the gold retaining screw is fastened 
to the abutment screw, clamping force on the implant is increased at the expense of 
decreasing the clamping force at the abutment screw-abutment interface by 50%. 
Maximum tensile stresses in the screw after preload were less than 55% of the 
yield stress. Cheong, et al.“ used FEA to predict that at a preload tension of 230N 
in the gold retaining screw shank, the clamping force at the abutment-abutment 
screw interface was first reduced to zero. With further tightening of the gold 
retaining screw, the rate of increase of stresses in it was faster than that of the 
abutment screw and it could be predicted that the gold retaining screw would fail 
first. Such failure by yielding was expected for a tensile load of around 400N 
applied to the gold cylinder. At this 400N tensile load, the clamping force at the 
fixture-abutment interface would be reduced to zero. This would affect the overall 
stability of the implant-prosthesis connection and eventually lead to component 
failure. 

Because preload application to the gold retaining screw reduces the clamping 
force at the abutment-fixture interface, it 1s recommended to find a balance preload 
between the gold retaining screw and abutment screw to make the whole implant- 
prosthesis connection more stable“. The current manufacturer recommendation for 
the Branemark system is to use tightening torques of 20N* cm for titanium 
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abutment screws and 10N * cm for gold retaining screws. 
3.3, 2.3 Washer 


The addition of a customized washer to dental implant screw joint systems may 
offer a very simple and inexpensive solution for the persistent problem of screw 
loosening. Versluis, et al.” studied the effect of a washer in a Branemark-type 
implant on the loosening conditions of the retaining screw using FEA. Their 
simulation indicated that a washer might significantly increase the axial tolerance of 
a screw against loosening up to 15 times more than a conventional system without 
washer. The authors state that this is accomplished by increasing the tolerance of 
the implant against deformation. 


3.3.3 Screw Fracture 


Factors that contribute to screw failure include the magnitude and direction of 
loading, the elastic modulus of the prosthesis and the rigidity of the abutment. 

Studying the IMZ implant system with FEA, Holmes, et al.” found that with 
increases in either load magnitude or load angle, stress concentrations in components 
of the implant system were generally increased. In another study, Holmes, et al.” 
also showed that in the IMZ implant, stress concentrations in bone and in 
components were much greater under a 30-degree load than under an equal vertical 
load. Greater deflection and stress concentrations within the coronal retaining screw 
were predicted with the use of a resin polyoxymethylene (POM) intramobile 
element (IME) than with a titantum element in the IMZ implant system. They also 
found that m FEA model stress transmission to bone was not reduced when the 
IME was modelled in POM rather than titanium. Maximum stress concentrations 
occurred in the fastening screw. 

Several authors***’ recommend high elastic modulus prostheses to avoid 
deflection of the prosthetic superstructure and stress concentrations in the retaining 
screw. Rigid abutment design is also needed to decrease the peak stresses in the 
screw and the deflection of the superstructure. Two related studies” described a 
FEA model of three different IMZ abutment designs: original threaded IME, 
Abutment Complete (ABC), and Intra-Mobile Connector (IMC). Progressive 
tightening of the retaining screw (preload) was simulated and the degree of screw 
tightening necessary to prevent opening of the crown-abutment interface in extreme 
loading (500 N occlusal load at 45 degrees) was determined individually for each 
system. A correlation was observed between the peak stresses in the screw and the 
deflection of the superstructure. Deflections and stress concentrations with the IMC 
were predicted to be in the same range as with the IME, but much greater than with 
the ABC. 


3.3.4 Summary 


The screw loosening problem is of concern especially when considering single tooth 
implant prostheses. The application of optimal preload has been the main means of 
preventing loosening. Moreover, a recent FEA study advocates the addition of a 
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washer as a simple and effective solution for the loosening problem. Stress 
concentration in fastening screws is influenced by load magnitude and direction. 
High rigid prostheses and abutments have been found to give more favorable stress 
distributions in screws. 


3.4 Multiple Implant Prostheses 


From a biomechanical viewpoint, there are three main classes of multiple wmnplant 
prostheses: (1) Implant-supported fixed prostheses (including cantilevered designs), 
(2) Implant-supported overdentures, and (3) Combined natural tooth and implant- 
supported prostheses. FEA studies for these prosthetic situations is usually more 
comp lex than for the single solitary implant. In most papers, three-dimensional FEA 
is considered to be necessary and two-dimensional FEA considered inadequate. 

Because multiple implants are splinted by the prosthesis framework, stress 
distribution is more complex compared to a single tooth implant situation. Loading 
at one point of the prosthesis will cause stress concentrations in all supporting 
implants to varying degrees. Also, the prosthesis can be loaded not by a single load 
but by multiple loads and in varying directions. In addition, the flexure of jaw bones, 
particularly the mandible, under functional loading conditions can cause stress tn the 
bone around the implants and may lead to bone resorption. Stress around the 
imp lant can be caused not only by local deformation of the bone due to movement 
of the implant and interface relative to the surrounding bone, but also by the 
complex deformation patterns of the mandible. 


3.4.1 Implant-supported Fixed Prostheses 


For implant-supported fixed prostheses, factors affecting bone-implant stress 
distribution and the ultimate success of prosthese, include implant inclination, 
implant number and position, prosthesis splinting scheme, occlusal surface and 
framework materials properties, and different cross-sectional beam shapes. 

Canay, et al.°° compared vertically orientated to angled implants and found that 
inclination of implants greatly influences stress concentrations around the implant- 
supported fixed prostheses. They found no measurable differences in stress values 
and contours when a horizontal load was applied to the vertical and angled implants. 
However, with vertical loading, compressive stress values were five times higher 
around the cervical region of the angled implant than around the same area in th 
vertical implant. | 

Many clinicians are of opinion that the selection of implant positions and the 
scheme of prosthesis splinting are critical for the longevity and stability of an 
implant prosthesis. Kregzde” reported that induced stresses in bone are sensitive to 
the scheme of prosthesis splinting and implant positions. He used three-dimensional 
FEA modelling of jaw bones, teeth and various implant numbers, positions and 
prosthesis designs to attempt optimisation of stress distribution to the implants. 
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The induced stresses on implants for different schemes of prosthesis splinting and 
different implant positions varied as much as 1,000% . 

The effect of different cross-sectional beam configuration for implant 
frameworks has also been investigated using FEA. Korioth and Johann” compared 
superstructures with different cross-sectional shapes and material properties during 
a simulated, complex biting task that modeled the deformation patterns of mandible 
during function. When they submitted their model to loads mimicking simultaneous 
bending and torsion of the mandibular corpus during a bilateral posterior bite, they 
found that predicted implant stresses varied significantly between implant sites for 
different superstructure shapes. The lowest principal stresses were obtained with a 
vertically orientated rectangular-shaped beam superstructure and contrary to 
expectations, the ideal “I-beam” superstructure cross section did not yield the 
lowest stresses. These authors concluded that implant abutment stresses were 
significantly affected by the cross-sectional shape of the prosthetic superstructure 
and by diverse mandibular loading conditions. 

Imp lant-sup ported fixed prostheses with cantilevers add additional factors that 
can influence stress distribution. These factors include cantilever length, cross- 
sectional beam shapes and recently, a system for additional support to the distal 
extension of the cantilever. Young, et al.” investigated a number of different cross- 
sectional beam shapes for cantilever fixed prostheses for initiation of permanent 
deformation at end loading. Straight and curved cantilever beams of 26 mm long were 
modeled in FEA. They found that a “L” shaped design was more rigid than other 
designs for a given mass, while an open “I” section framework offered good 
possibilities particularly when used as curved shapes. “L” shaped cobalt-chromium 
or stainless steel frameworks of 26 mm cantilever span underwent permanent 
deformation at end loadings between 130 and 140 N depending on section curvature. 
They caution that a good framework design is critical to avoid failures since it is 
known that biting loads can exceed these values. 

Different material properties affect stress distribution in different ways. Korioth 
and Johhann” showed that an increase in elastic modulus of prosthetic materials 
does not necessarily lead to a decrease in stresses on all existing implant abutments. 
Less rigid superstructures seemed to increase implant abutment stresses overall as 
well as to decrease tensile stresses on the most anterior implant abutments in the 
modeled complex biting task. 

Using a six implant-supported mandibular complete arch fixed prosthesis three- 
dimensional FEA model, Sertgdz investigated the effect of different occlusal surface 
materials (resin, resin composite, and porcelain) and different framework materials 
(gold, silver-palladium, cobalt-chromium, and titanium alloys) on stress distribution 
in the fixed prosthesis and surrounding bone. He demonstrated that using a 
prosthesis superstructure material with lower elastic modulus did not lead to 
substantial differences in stress patterns or levels in the cortical and cancellous bone 
surrounding the implants. For a single loading condition investigated, the optimal 
combination of materials was cobalt-chromium for the framework and porcelain for 
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the occlusal surface. 

Sertgdz and Guvener™, with a three-dimensional FEA model of a bilateral distal 
cantilever fixed prosthesis supported by six implants in the mandible, predicted that 
maximum stresses occurred at the most distal bone-implant interface on the loaded 
side and this significantly increased with increase in cantilever length. Nevertheless, 
they found no significant change in stress levels associated with implant length 
variation. However, Lindquist, et al.” in a 15-year longitudinal clinical follow-up 
study, reported that bone at the distal implants of cantilevered mandibular imp lant- 
supported prostheses remained very stable and conversely, more bone loss was 
observed around the anterior implants. However, this may be caused by a multitude 
of clinical factors. It was concluded that occlusal loading factors such as maximal 
bite force, tooth clenching and cantilever length were of minor importance to bone 
loss in their study population. This suggests that extrapolation of FEA studies to 
clinical situations should be approached with caution. 

New systems for additional support to the distal extensions of cantilevered 
prostheses have been suggested. IL system uses a short implant and a special ball 
type attachment to support the distal extension of cantilevered prostheses. 
Lewinstein, et al.” compared this new support system to a conventional cantilever 
prosthesis using two-dimensional FEA. This system dramatically lowered the 
stresses in bone, cantilever, and implants potentially reducing failures within the 
implants, prostheses, and surrounding bone. Moreover, employment of a relatively 
long-span prosthetic extension in the posterior region of the jaw would be possible. 

In summary, stress distribution in implant-supported fixed prostheses have been 
shown by FEA to be influenced by the factors of implant inclination, implant 
number and positions, prosthetic splinting scheme, superstructure material 
properties and beam design. 


3.4.2 Implant-supported Overdentures 


The use of implant-supported overdentures is viewed as a cost-effective treatment 
modality. Some clinicians are of the opinion that the designed stress-breaking 
features of overdenture attachments confer more favourable biomechanical 
characteristics compared to implant-supported fixed prostheses. Implant-supported 
overdenture attachment systems include bar-clips, balls, O-rings and magnets. For 
bar-clip attachment systems, some biomechanical factors identified are the number 
of implants, bar length, stiffener height, and material properties. 

Meijer, et al.” set up a three-dimensional model of a human mandible with two 
endosseous implants in the interforaminal region and compared stress distribution 
when the two implants were connected by a bar or remained solitary. The most 
extreme principal stress was found with oblique bite loads whereas vertical bite 
loads resulted in the lowest stress. The most extreme principal stresses in bone were 
always located around the necks of the implants. No significant differences in stress 
distribution were predicted with the highest maximum and lowest minimum 
principal stresses being 7.4-and-16.2 MPa in the model without the bar and 6.5-and- 
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16.5 MPa with the bar. It was also found that a bar placed anteriorly of the 
interconnecting line between the two implants caused extremely large compressive 
and tensile stress concentrations in the bone around the implants. Therefore, in 
those cases, it is advised not to connect the implants or, if a bar-clip attachment is 
preferred, to place additional implants in the frontal region’. In a further study, 
Meijer, et al. used the same model to study a four implant system with the 
implants either connected by a bar or remaining solitary. Their results showed that 
with uniform loading, there were more or less equal extreme principal stresses 
around the central and lateral implants. With non-uniformly loading on the 
superstructures, the implant nearest to the load showed the highest stress 
concentration and with connected implants there was a reduction in the magnitude of 
the extreme principal stresses compared to solitary implants. 

FEA modelling of a two implants round bar” and Hader bar system” as well as 
a four implant Hader bar system® found span length and stiffener height to be more 
critical factors in the adequacy of the overall design as compared with changing 
material properties in the range of alloy stiffness tested. 

Overdentures supported by two implants ball system have achieved better 
stress distribution in bone compared to a two implants bar system. Menicucci, et 
al.” used three-dimensional FEA to evaluate transmission of masticatory load in 
mandibular imp lant-retained overdentures. Overdentures retained either by two ball 
attachments or by two clips on a bar connecting two implants were compared. For 
the ball attachment system, a 35 N load on the first mandibular molar of the 
overdenture induced a greater reaction force on the distal edentulous ridge mucosa of 
the nonworking side compared to the bar-clip attachment. However, when pert 
implant bone stress was considered, this was greater with the bar-clip attachment 
than with the ball attachment. 

In summary, FEA has been used to investigate the stress distribution obtained 
when implants were left solitary, used with ball attachments or connected by bars 
for clip retention in various configurations and designs. Not all studies modelled the 
overdenture over the implants and bar superstructure. Bar design factors like 
stiffener height and span length were found to significantly affect stress distribution 
whereas the influence of various material moduli was comparatively less significant. 


3.4.3 Combined Natural Tooth and Implant-supported Prostheses 


Combining natural teeth and implants to support fixed prostheses has been 
advocated by certain quarters in implant dentistry. Controversy exists as to the 
advisability of this design philosophy from a biomechanical as well as a clinical 
perspective. A significant clinical consideration in the restoration of partial 
edentulism with implant and tooth-supported prostheses is whether implants and 
natural teeth abutments should be splinted, and if so, in what manner. There is a 
differential deflection between the viscoelastic intrusion of a natural tooth in its 
periodontal ligament and the almost negligible elastic deformation of an 
osseointegrated implant. This difference may induce a fulcrumrlike effect and 
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possibly overstress the implant or surrounding bone. The biomechanical factors that 
can influence the stress distribution include abutment design, implant material 
properties, effect of resilient elements, connector design (precision or semi-precision 
attachments), and degree of splinting implants to natural tooth abutments. 

For the implant connected with a natural tooth situation, van Rossen, et al.” 
concluded that a more uniform stress was obtained around implants with stress- 
absorbing elements of low elastic modulus. It was also concluded that the bone 
surrounding the natural tooth showed a decrease in peak stresses in such a situation. 

Charkawi, et al.” studied the use of a resilient layer material under the 
superstructure of the implant in a connected tooth and implant supported 
prosthesis model. Their FEA results proposed that this new modification could 
mimic the structural natural tooth unit by allowing movement of the superstructure 
without movement of the implant when loaded. 

However, Misch and Ismail* conducted a three-dimensional FEA comparing 
models representing a natural tooth and an integrated implant connected by rigid and 
nomrigid connectors. Based on similarities in stress contour patterns and the stress 
values generated in both models, they concluded that advocating a non-rigid 
connection because of a biomechanical advantage may be erroneous. Melo, et al.” 
also investigated tooth and implant-supported prostheses in free-end partially 
edentulous cases. Their two-dimensional FEA predicted that lowest levels of 
stresses in bone occurred when the prosthesis was not connected to a natural 
abutment tooth but was supported instead by two free-standing implant abutments. 
Normrigid attachments, when incorporated into a prosthesis, did not significantly 
reduce the level of stresses in bone. A recent comprehensive review of both clinical 
and laboratory studies concluded that the issue of connecting natural teeth to 
implants with rigid or non-rigid connectors still remains unresolved*’. 


3.5 Conclusions 


FEA has been used extensively in the prediction of biomechanical performance of 
dental implant systems. This chapter has reviewed the use of FEA in relation to the 
bone-implant interface, the implant-prosthesis connection, and multiple implant 
prostheses. 

Load transmission and resultant stress distribution at the bone-implant interface 
has been the subject of FEA studies. Factors that influence load transfer at the 
bone-imp lant interface include the type of loading, implant and prosthesis material 
properties, the length and diameter of implants, implant shape, the structure of 
implant surface, nature of bone-implant interface, and the quality and quantity of 
surrounding bone. Of these biomechanical factors, implant length, diameter, and 
shape are easily modified in the implant design. Cortical and cancellous bone quality 
and quantity needs to be assessed clinically and can influence implant selection. 

Stress distribution in the implant-prosthesis connection has been examined by 
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FEA studies because of the incidence of clinical problems such as gold screw 
failures, abutment screw failures, and implant fracture. Design changes to avoid or 
reduce these prosthetic failures by improving the stress distribution of implant 
components have been suggested. 

When applied to multiple implant prosthesis design, FEA has suggested 
improved biomechanical situations when factors such as implant inclination, implant 
positions, prosthetic material properties, superstructure beam design, cantilever 
lengths, bar system, bar span length and stiffener height, and overdenture attachment 
type were optimised. For combined natural tooth and implant supported 
prostheses, FEA studies were inconclusive whether to use a rigid or resilient 
implant systems. 

FEA is an effective computational tool that has been applied from the 
engineering arena to dental implant biomechanics. Many design feature optimisations 
have been predicted and will be applied to potential new implant systems in the 
future. 
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1.1 Introduction 


The finite element method may be applied to all kinds of materials in many kinds of 
situations: solids, fluids, gases, and combinations thereof; static or dynamic, and, 
elastic, inelastic, or plastic behaviour. In this book, however, we shall restrict the 
treatment to the deformation and stress analysis of solids, with particular reference 
to dental implants. 


1.2 Historical Development 


Deformation and stress analysis involves the formulation of force-displacement 
relationships. These have been used in increasingly sophisticated forms from the 
1660s, when Robert Hooke came out with his Law of the Proportionality of Force 
and Displacement. 

The nineteenth and twentieth centuries saw a lot of applications of the force 
displacement relationships for the analysis and design of large and complex 
structures, by manual methods using logarithmic tables, slide rules, and in due 
course, manually and electrically operated calculators. 

Particular mention must be made of the contributions of the following scientists, 
relevant to modern structural analysis: 

1857: Clapeyron Theorem of Three Moments 

1864: Maxwell Law of Reciprocal Deflections 

1873: Castigliano Theorem of Least Work 

1914: Bendixen Slope-deflection Method 

References for these works and others to follow are given at the end of the 
chap ter. 

These and other early methods and applications to articulated (stick-type) 
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structures were based on formulas developed from structural mechanics principles, 
applied to straight, prismatic members such as axial force bars, beams, torsion rods, 
etc. 

All these techniques yielded simultaneous equations relating components of 
forces and displacements at the joints of the structure. The number of simultaneous 
equations that could be solved by hand (between 10 and 15) set a practical limit to 
the size of the structure that could be analysed. 

To avoid the direct solution of too many simultaneous equations, successive 
approximation methods were developed. Among them should be cited the following: 

1932: Cross Moment Distribution Method 

1940: von Karman and Biot Finite Difference Methods for Field Problems 

1942: Newmark Finite Difference Methods for Structural Problems 

1946: Southweil Relaxation Methods for Field Problems 

These expanded the size limitations outwards by many orders of magnitude, 
enabling large complex articulated as well as plate-type structures to be analysed 
and designed. 

The appearance of commercial digital computers in the 1940s revolutionised 
structural analysis. The simultaneous equations were not an obstacle any more. 
Solutions became even more efficient when the data and processing were organised 
in matrix form. Thus was matrix analysis of structures born. 

It was the aeronautical industry that exploited this new tool to best advantage, 
but structural designers were quick to follow their lead. By the 1960s, not only 
could better and bigger aircraft be manufactured, but large bridges and buildings of 
unconventional design could be built. 

This also resulted in the computerised revival of the somewhat abandoned earlier 
methods of consistent deformation and slope deflection. Not only could much larger 
problems be handled, but also effects formerly neglected as secondary (out of 
computational necessity) could be included. Pioneers in matrix computer analysis 
were: 

1958: Argyris- Matrix Force or Flexibility Method 

1959: Morice-Matrix Displacement or Stiffness Method 

From matrix analysis of articulated structures to finite element analysis of 
continuous systems, it was a big leap, inspired and spurred on by the digital 
computer. However, it was not as if the entire idea was new. 

Actually, the history of the Finite Element Method is the history of 
discretisation, the technique of dividing up a continuous region into a number of 
simple shapes. The progress from conceptualisation and formalisation, to 
implementation and application, may be summarised as follows: 

1774: Concepts of Discretisation of Continua (Euler) 

1864: Framework Analysis (Maxwell) 

1875: Virtual Work Methods for Force-displacement Relationships (Castigliano) 

1906: Lattice Analogy for Stress Analysis (Wieghardt) 

1915: Stiffness Formulation of Framework Analysis (Maney) 
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1915: Series Solution for Rods and Plates (Galerkin) 

1932: Moment Distribution Method for Frames (Hardy Cross) 

1940: Relaxation (Finite Difference) Methods (Southwell) 

1941: Framework Method for Elasticity Problems (Hrenikoff) 

1942: Finite Difference Methods for Beams and Columns (Newmark) 

1943: Concept of Discretisation of Continua with Triangular Elements (Courant) 

1943: Lattice Analogy for Plane Stress Problems (McHenry) 

1953: Computerisation of Matrix Structural Analysis (Levy) 

1954: Matrix Formulation of Structural Problems (Argyris) 

1956: Triangular Element for Fintte Element Plane Stress Analysis (Turner, et al.) 

1960: Computerisation of Finite Element Method (Clough) 

1964: Matrix Methods of Structural Analysis (Livesley) 

1963: Mathematical Formalisation of the Finite Element Method (Melosh) 

1965: Plane Stress and Strain, and Axisymmetric Finite Element Program 

(Wilson) 
1967: Finite Element Method in Structural and Continuum Mechanics 
(Zienkiewicz) 

1972: Finite Element Applications to Nonlinear Problems (Oden) 

Old theories of solid continua were reexamined. Up to the 1950s, only 
continuous uniform regions of some regular shape such as square and circular plates 
or prisms could be analysed with closed form solutions. Some extensions were made 
by conformal mapping techniques, Series and finite difference solutions were 
developed for certain broader class of problems. But all these remained in the 
domain of academic pursuit of theoretical advancement, with few general 
applications and limited practical use. 

Again, it was the aircraft industry that pioneered the idea of analysing a region 
as the assemblage of a number of triangular elements. The force-displacement 
relationships for each element were formulated on the basis of assumed 
displacement functions. The governing equations resulted after approximately 
assembly modelled the behaviour of the entire region. Once the equations were 
formulated, further solution followed the same steps as the matrix structural 
analy sis. 

The idea worked, and very efficiently with computers. It was also confirmed 
that the finer the division, the better the results. Now the aircraft designers could 
consider not only the airframe, but the fuselage that covered it and the bulkheads 
that stiffened it, as a single system of stress bearing components, resisting applied 
forces as an integrated unit. 

This technique came to be called the “Finite Element Method” (“FEM” for 
short), both because a region could be only broken up into a finite number of 
elements, and because many of the ideas were extrapolated from an infinitesimal 
element of the theory to a finite sized element of practical dimensions. 

Clough and his associates brought this new technique into the civil engineering 
profession, and soon engineers used it for better bridges and stronger shells. 
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Mechanical engineers exploited it for understanding component behaviour and 
designing new devices. 

Computer programs were developed all over the Western world and Japan. The 
first widely accepted program was “SAP” (for Structural Analysis Package) by ELL. 
Wilson, which got him a Ph.D. from the University of California, USA. Most 
programs were in FORTRAN, the only suitable language at the time. Soon there was 
a veritable explosion in programs, and today, there are scores of packages in recent 
languages which are menu-driven and automated to the extent that with minimal 
(self-)training, anybody can do a finite element analysis for better or for worse! 

Purists viewed the early applications with considerable reservation, pointing out 
the lack of mathematical rigour behind the technique. Appropriate responses were 
not slow in coming. Melosh and others soon connected the assumptions behind the 
formulation of the element with the already prevalent classical methods of 
interpolation functions. 

Argyris in Europe, Zienkiewicz in UK, and Clough, Wilson, Oden, and numerous 
others in USA, pushed the frontiers of finite element knowledge and applications 
fast and wide. Between the 1950s and 1970s, applications of the finite element 
method grew enormously in variety and size, supported or triggered by fantastic 
developments in digital computer technology. In the last two decades, new 
developments have not been so many, but practical applications have become wider, 
easier, and more sophisticated. 

Early users, the author included, considered hundred elements as a boon. A 
decade later, third generation computers enabled analysts to routinely use thousands 
of elements. By the 1970s, capacity and speed had increased ten times further. 
Nothing seemed to be beyond reach of finite element analysis whether it be a 
nuclear reactor (Fig. 1.1(a),(b),(c)), or a tooth (Fig. 1.1(d)), both of which the author 
has analysed. 
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Fig. 1.1 (a) Test Model of Prestressed Concrete Nuclear Reactor; (b) One-twelfth Symmetry 
Segment for Analysis; (c) 3-D Finite Element Idealisation of the Analysis Segment; (d) 3-D 
Finite Element Idealisation of a Tooth 


Now, computer packages which once demanded a mainframe have come to the 
desk top, and been loaded with powerful program graphics user interfaces, and 
interactive, online modelling and solutions. 
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It was just a small imaginative step to extend the applications beyond linear 
structural analysis, to non-linear and plastic behaviour, to fluids and gases, to 
dynamics and stability, to thermal and other field problems, because all of them 
involved the same kind of differential equations, differing only in parameters and 
properties, while the overall formulation, assembly, and solution techniques 
remained the same. 

The references of historical importance, given at the end of the chapter, are 
merely representative, often the earliest in a series of many publications on a topic 
by the same or other authors. More detailed coverage of the history and further 
references may be found in the works by Cook, Desai, Gallagher, Huebner, Oden, 
Przemieniecki, and Zienkiewicz. Readers can referr to these resources for additional 
information on any of the topics discussed by the author in the following chapters. 

Today, there is almost no field of engineering, no subject where any aspect of 
mechanics is involved, in which the finite element method has not made and is not 
continuing to make significant contributions to knowledge, leading to unprecedented 
advances in state of the art and its ultimate usefulness to mankind including 
contributions to dentistry. 


1.3 Definitions and Terminology 


The basic procedure for matrix analysis depends on the determination of 
relationships between the “Actions”, namely forces, moments, torques, etc. acting 
on the body, and the corresponding “Displacements”, namely deflections, rotations, 
twists, etc. of the body. 

A “structure” is conventionally taken to consist of an assembly of straight 
“members” (as in trusses, frames, etc.) or curved lines whose shape can be 
mathematically evaluated, which are connected, supported, and loaded at their ends, 
called “joints”. Fig.1.2(a) shows a two-storey structure consisting of frames in the 
vertical plane, grids in the horizontal plane, and trusses for the entrance canopy. 

A “system” conventionally consists of a continuous membrane, plate, shell, or 
solid, single or in combination, each divided for analysis purposes into a finite 
number of “elements”, connected, supported, and loaded at their vertices and other 
specified locations on edges or inside, called “nodes”. Systems may include 
structures as well. 

Fig. 1.2(b) shows a machine part system consisting of a solid, thin-walled shell, 
and a projecting plate. The suggested divisions are shown tn lines of a lighter shade. 
Generally, the curved boundaries will be modelled as straight lines. The circular pipe 
in this case will be simulated as a hexagonal tube. 

The principal difference between a structure and a system is this: The 
articulated structure is automatically, naturally, divided into straight (and certain 
regularly curved) members such as the truss member AB in Fig. 1.2(a), whose 
behaviour is well known and can be formulated theoretically. On the other hand, the 
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Fig. 1.2 (a) Tworstorey Articulated Structure; (b) Machine Part System, Continuum 


continuous system has no such theoretical basis and has to be divided into pieces of 
simple shape, such as the triangle IJK in Fig. 1.2(b), whose behaviour must be 
formulated by special methods. 

Most real life facilities involve a combination of both types described above. For 
instance, a building has flat platetype walls and floors; a machine may sit atop 
columns and beams. In practice, “member” and “joint” usually apply to a structure, 
while “element” and “node” apply to a system in particular, and to a structure also 
in general. 

Each node or joint can have a number of independent action (force or moment) or 
displacement (deflection or rotation) components called “Degrees Of Freedom” 
(DOF) along a certain direction corresponding to coordinate axies. 

A plane truss member such as AB in Fig. 1.2(a) shown enlarged in Fig. 1. 3{a) 
has two DOF at each joint. Hence the member has a total of (2X 2) or 4 DOF. 


(a) 6) B (b) o 
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Fig. 1.3 (a) A Truss Member AB; (b) A Triangular Finite Element IJK 


A triangular membrane element such as IJK in Fig, 1.2(b) shown enlarged in Fig, 
1.3(b) has two DOF at each node. Hence the element has a total of (3 X 2) or 6 DOF. 

Different types of members and elements have different numbers of DOF. For 
instance, a 3D frame member has two joints and six DOF (3 forces or displacements 
and 3 moments or rotations) per joint and 12 DOF in total. A solid “brick” element 
has eight nodes and three DOF (3 forces or displacements) per node and 24 DOF in 
total. 

Additionally, in the case of finite elements, joint the same type of element may 
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even have different number of nodes in “transition” elements. 


1.4 Flexibility Approach 


Fig. 1.4 shows a truss member with actions and corresponding displacements along 
the two DOF at each end. The sets of four actions and displacements can be 
represented vectorially or in terms of x, y components, as follows: 


{At = {A,, Ao, Aj, As} = {X;, Yi, Xj, Y;}, the “Action Vector” 
{D}= {D,, D,, D3, Dy} = {ui, vi, uj, vj}, the “Displacement Vector” 


The displacement D at every DOF (say J) is a function of the actions A,, A:, at 
all connected DOF. Within the elastic limit, D; is a linear combination of the effects 
of all actions. 

Thus, their relationship may be written as: 


D,= fi, Ait fiz Ao+ fis Ast fia Ag (1.1) 


in which f; stands for the displacement at DOF I 
due to a unit action at DOF J, and is known as the 
“Flexibility Coefficient”. f) 4) 

Three more such equations may be written for ~ ‘__- D, 
D,, D;, and D,. The four equations may be 
represented in matrix form as: 


{D}=LF] + {A} (1.2) 


4X1 4X4 4X1 


Fig. 1.4 Displaced Truss 
Member 


in which, the [F] matrix of flexibility coefficients 
is known as “Flexibility Matrix”. 

The flexibility coefficients for prismatic bars can be determined from basic 
theoretical principles of strength of materials and theory of structures. 

The flexibility approach was quite popular as the “Force Method” for manual 
analy sis, the “Method of Consistent Deformation” being a typical application. With 
the advent of computers, it was found that this approach was not convenient to 
formulate or solve large and complex problems. Hence, the flexibility approach was 
not pursued further for practical applications. 


1.5 Stiffenss Formulation 


1. 5.1 Stiffness Matrix 


An alternative formulation, an exact opposite—in fact the inverse— of the flexibility 
approach, called “stiffness approach” or “displacement approach” was also in use 
for manual solutions. The “Slope Deflection Method” for continuous beams and the 
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“Moment Distribution Method” for beams and frames were very popular. 

This approach was very convenient for computerisation and became the 
preferred method for computer solutions, especially for finite element analysis. 

In general, the displacement along every DOF needs an action along that DOF 
and reactions at all the other connected DOFs for equilibrium. For elastic behaviour, 
the function 1s a linear combination of all the displacement effects. 

Thus, the action-displacement relationships of the truss member in Fig, 14 is 
written as: 


A. = ky D,+ ky» D.+ ky; D3+ ki Dy 
Ao= ky, D, + ky Dy + ky; D3 + Ky Dy 
A3= ky, D, + ky. D. + ky; D3 + ky, Dg 
Ag= Kg Dy + Ky D+ kag Dy + Kay Dy (1.3) 


in which k; stands for the action at DOF I due to a unit displacement at DOF J 
(with all other displacements set to zero) and is known as the “ Stiffness 
Coefficient”. 

- The four Eg.(1.3) may be represented in matrix form as: 


{A}=Lk] {D} (1.4) 


AX] 4X4 4X1 


in which the [k] matrix of stiffness coefficients is known as “Stiffness Matrix”. 
The stiffness coefficients for prismatic bars can be determined from basic 
theoretical principles of strength of materials and theory of structures. 
For instance, consider the truss member AB, of length L and cross-sectional area 
A, from a material with Young's Modulus of elasticity E, inclined at an angle @ with 
the horizontal, subjected to a unit displacement along DOF 1, as shown in Fig. 1.5 


(a). 

A,=A ,sin@ 
(b) 
A,=A,cos@ 


A,=A,cos@ 


al D =! 


=kcos0” A,= A,sin@ 


Fig. 1.5 (a) Unit Global Displacement; (b) Action Components 


The unit horizontal displacement D, resolves into an axial displacement D,= 1 * 
cosé = cos@ and a transverse displacement D,’= 1 ° sin@ = sind. 

As the truss member ends are pinned, only the axial displacement D, needs a 
force A,= kD,, or kcos@, k being the stiffness of the axial force bar, namely (EA,/L). 

As shown in Fig. 1.5(b), this axial force A, may now be resolved into: 
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A,= A,cos@ =kcos’@ and A,= A,sin? = kcos@sin@ 


To keep the bar in equilibrium, equal and opposite reactions must be developed 
at the end B, giving: 


A;= — A,;=—kcos’*@ and A= — A.= — kcos@sind 


These four actions corresponding to a unit displacement along DOF | are defined 
by the four stiffness coefficients kj (i= 1, 2, 3, 4) that constitute column 1 of the 
stiffness matrix [K] in Eq. (1.4), as shown in bold type in Eq. (1.5) below. The other 
three columns can similarly be determined by the application of unit displacements 
along DOF 2, 3, and 4 in turn. 

The resulting stiffness matrix is as follows: 


Ay cos’ cosdsing — cos’ — cos@sind } | D, 

A, |_ EA, | cos@sin@ sin’ —cosésind —sin’é D; (1.5) 
A; L |—cos’6 —cosésin? cos’é cos@sing D, 
A, —cosésing? —sin’é sindcos? sin’ 6 D, 


Stiffness matrices can be developed for other straight prismatic members such as 
beams and torsion bars from similar principles. 

However, the situation is quite different when it comes to finite elements. 

The triangular plane element IJK in Fig. 1.6, 
under the action of six force components along the 
6 DOF, ts represented as: 


{A}= {A1, A, As, Aa, As, As} 
= {X,, Y:, X;, Yj, Xx, Yx} also. 


It is displaced to the configuration I J’K’, with 
the deflection components: 


{D}= {D,, D,, D;, Da, Ds, Ds} 
= fui Vis Uj, Vj, Uk, Vi } also. 


Fig. 1.6 Action and Displacement 


Relationships of actions and displacements at 
Components 


the DOF of this element are of the same kind as 
Eqs. (1.3) and (1.4) for the truss member, with the 
difference that for the element, the {A} and {D} vectors are (6X 1) and stiffness 
matrix [K] is (6X 6). 

However, it is unlike Eq. (1.5) that no theoretical method to determine the 
stiffness coefficients for a general triangle or any other shape exisxs. Other special 
techniques must be resorted to, as will be discussed in subsequent chapters. 


1.5, 2 Characteristics of Stiffness Matrix 


The characteristics of the member or element stiffness matrix, most of which may be 
deduced from Eq. (1.5), are listed below as common to all element stiffness matrices. 
(1) The stiffness matrix is square, logically from the fact that there are as many 
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displacement DOF as action DOF. 

(2) The stiffness matrix is symmetric. This derives from the principle of 
conservation of energy, commonly developed as Maxwell's Law of Reciprocal 
Deflections for structural members, which states that the displacement at A due 
to a unit action at B is equal to the displacement at B due to a unit action at A. 

(3) The matrix is “positive definite”, that is the diagonal terms are positive and 
(generally) dominant. It simply reflects the fact that a point at which an action 
is applied, moves along the direction of the action, not opposite to it. 

(4) Each column of the matrix, representing all the actions acting on the element due 
to the displacement at one DOF, must satisfy statics. If they are all forces, then 
the alternate terms in a row or column (representing horizontal and vertical 
components separately) must add up to zero separately. Moments of all the 
forces about any point must be zero. 

(5) The determinant of the stiffness matrix will be found to be zero. This implies 
that the matrix is singular, and cannot be inverted. In effect this means that the 
displacements due to any action on the member will be infinite, that is, not 
capable of being determined. 

The explanation for this apparent anomaly is very simple: The stiffness matrix 
is just a property of the element. The element can accept and resist a load only 
when it is supported against rigid-body movement, as along three DOF in 2D. 
Without such minimal support, even the smallest load along any of its DOF will 
simply blow the element away to “infinity” as in space! 


1.5.3 Equivalent Loads 


Loads are often applied between the joints of a member, such as a distributed load 
on a beam. As matrix analysis deals with loads and displacements at only the joints, 
the member loads must be replaced by “Equivalent Loads” at the joints. These 
actions are also called “Consistent Nodal Loads”. 

In this case, for straight prismatic members, classical theories provide values for 
equivalent loads. 

For instance, the beam of span L loaded with q per unit length shown in Fig, 1.7(a) can 
be replaced with the two forces and two moments shown in Fig. 1.7(b), on the basis 
that both of them produce the same end rotations @, and satisfy statics. 


(a) q 


Fig. 1.7 (a) Simply Supported Beam with Uniform Load; (b) Equivalent End Actions 


Situations in regard to finite elements are not as simple as this and will need 
special treatment. 
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1.5.4 System Stiffness Equations 


For a system with n DOF, the governing equation for the I-th DOF of the 
assemblage of members or elements is obtained by combining the governing 
equations for the same DOF from the individual pieces, in the form: 


Aj= ky D+ kp, D,+°°°+ kD, 
or, for all the n DOF, in matrix form, similar to the element Eq.(1.4): 
{A,} =[K,}{D,} (1.6) 


where, the action vector {A,} includes the effects of internal loads and is (nX 1) in 
size; the displacement vector {D,} is (nX 1) in size; and the stiffness matrix [K,] is 
(nX n) in size. 

Since the system stiffness matrix is the superposition of the element stiffness 
matrices, all the characteristics of the element stiffness matrix listed in Section 1.4.2 
can carry over into the system stiffness matrix. 


1. 6 Solution Methodology 


Typically, the data for a problem in structural and continuum analy sis consists of: 
(1) Geometry, namely location of nodes; 

(2) Topology, namely the nodes by which various elements are connected; 

(3) Relevant material properties; 

(4) Locations of supports, and their movements, if any; 

(5) Locations and magnitudes of loads. 

To solve a finite element problem, first the region is divided into sufficient 
numbers of elements of a suitable type, to reflect the geometry and any special 
features, with nodes located at supports and concentrated loads. 

From the data on geometry, topology, and material properties, a stiffness matrix 
for each element is determined. 

The element stiffness matrices are assembled to form the system stiffness 
matrix, as will be explained in the next chapter. 

The supported and loaded DOF in data items (4) and (5) form a complementary 
set: Those DOF which are constrained can develop reactions and hence must not be 
loaded. Differently, any action component that happens to be applied on a 
supported DOF will directly pass on to the support and hence must be excluded 
from the main analysis. Likewise, those DOF that are loaded (and unconstrained) are 
free to displace. 


1. 6. 1 Manual Solution 


Let a be the number of known actions and unknown displacements and c the number 
of unknown reactions and known displacement constraints. Note that total DOF n 


12 Application of the Finite Element Method in Implant Dentistry 


=artc. 

Thus the action and displacement vectors can be partitioned in a mutually 
exclusive manner as {A,| A.} and {D,| D.}. Note that the action vectors are actually 
modified action vectors incorporating the effects of any member or element loads. 

Subsequently, the stiffness equations can be rearranged to match such 
partitioning and the stiffness matrix can also be partitioned into four parts. 

The partitioned sub-vectors and sub-matrices and their sizes, are as follows: 


A, K. ‘Kz D, 

(aXl){_ | @Xa) (@XQ] J (aX) Gia) 
Ae Re a Kee Dz 

(cX1) cXa)' (cXe (cX}) 


This can be separated into two matrix equations as follows: 


{A.}= [Ka ]{D.}+[K. ]{D.} (1.8a) 
and {Ac }= [Ka ]{D.}+[Ke ]{D-} (1.8b) 


In Eq. (1.8a), all the terms except {D,} are known, and {D,}can be computed 
from: 


{Aa }-([Ke ]{D.}= [Ka ]{D;} 


Le. {D.}= [Ka] ({A.}([Ke ]{D-}) (1.9a) 
Designating the term (-[K,.]{D.}) as {Aa}, Eq. (1.9a) may be written as: 
{D.}= [Ka] '({Aa}+ {Aa}) (1.9b) 


If all known (support) displacements {D,} are zero, {Az} is zero. Eqs. (1.9) 
simplify to: 


{D.}=[Ka]' {Ao} (1.10a) 


Now, with all the displacements known, the unknown reactions {A,} may be 
found from Eq. (1.8b). 

Further, if all the known (support) displacements {D,} are zero, then {Aj} is 
zero, and Eq. (1.8b) simplifies to: 


{A.}= [Ka ]{D,} (1.10b) 


1.6.2 Computer Solution 


While the partitioning procedure described may be suitable for manual solution of 
small problems, obviously the reordering of the terms in, and partitioning of, large 
vectors and matrices in computer solutions will be inefficient and time consuming. 

Instead, the equation set is retained as assembled. The action vector {A,} is still 
modified to {A;} for the c known displacements according to Eq. (1.9a), but 
without rearranging the equations. The stiffness coefficients k; in the columns 
corresponding to the known displacements are set to zero as their effect has already 
been incorporated into the action vector. 


1 Finite Element Method 13 


As those a equations are the only ones needed for solution of the unknown 
displacements, the other c equations corresponding to Eq. (1.8b) involving the 
known displacements and unknown reactions are replaced by dummy equations in 
the form of {A; }=[K.]{D.}, in which each A; is set equal to the known value of 
kK'D,, so that when D; is solved for, it will return the known value of D,. 

The equations thus modified, still in their original order, can be written as: 


{A’}=[K']{D"} (1.11) 


Because of the incorporation of the support conditions, the modified [K’] matrix 
is not singular any more. Hence the displacements may be found as: 


{D'}=[K'}' {A} (1.12) 


Then all the displacements can be computed, including the known displacements 
at their original input values. Once the displacements are known, the unknown 
reactions can be computed from the full set of original equations, Eq. (1.6). 

Needless to say, the computer ts also heavily involved in the automation of the 
input preparation and output evaluation of such large scale analyses. 


1.6.3 Support Displacements 


It may be noted that Eqs. (1.9) are in general form wherein some or all the known 
displacements may be non-zero, implying support settlement or yielding. 

The minimum supports that a system must be provided before analysis can 
proceed is strong enough to prevent rigid body displacement. For instance, in a 2D 
plane region, three non-collinear DOF must be supported to prevent rigid body 
deflection and rotation. 

In such a minimally supported system, any support displacement will only 
cause a change in the position of the body and will not introduce deformations. 
Hence no internal actions will be developed due to support displacement. The 
reactions at the supports can be found from statics, and the only internal actions 
will be due to external applied loading, if any. 

However, as is more common, if the system is supported at more than the 
minimum required number of DOF, then it becomes “statically indeterminate”. Any 
support displacement will introduce internal deformations and actions, even without 
external applied loading. Eqs.(1.7), (1.8), and (1.9) will take care of all these effects. 

The notation {A;}=[K,]|{D.} introduced in Eqs. (1.9), may now be interpreted 
as the “Equivalent Load” vector to account for support displacements. 


1.6.4 Alternate Loadings 


Note that Eqs. (1.9) involve the applied loading conditions {A,} on the right hand 
side only. Hence if we need the results for different applied loadings, we can simply 
save the [K.]' matrix and the {A,} vector, and carry out the matrix multiplication 
in Eqs. (1.9) for the new applied load vector {A,}. 

From the displacements {D,} due to the new loading, the corresponding 
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reactions, {A,} may be found from Eq. (1.8b) or Eq. (1.10b). 

This facility is of great use when different loadings have to be applied to the 
same object in the same support conditions, as is very often the case. 

Thus, if the results can be computed for a few basic independent loadings, then 
for various combinations of these loadings, the stiffness equations do not have to be 
solved afresh. Only the already computed results of the basic loadings need to be 
combined in the same proportions as the loads in the combinations. 

As should be evident, the solution methodology described is standard practice 
and quite routine with today’s computers. The art and science of finite elements 
thus revolves entirely around the determination of the stiffness matrix and the 
equivalent load vector. 


1.7 Advantages and Disadvantages of FEM 


The advantages of FEM, some already touched upon, may be summarised as 

follows: 

(1) Any domain with curved boundaries, heterogeneous material properties, irregular 
support constraints, and varying loading conditions, may be sub-divided into a 
suitable number of fmite elements, appropriate material and behaviour properties 
may be ascribed to them, and the resulting governing equations may be solved 
quickly and accurately by computers. 

(2) It is equally applicable to statics and dynamics; solids, fluids, and gases and 
combinations; linear and nonlinear; elastic, inelastic, or plastic; special effects 
such as crack propagation; events and processes such as bolt pretensioning, etc. 

(3) The massive amounts of data itself can be efficiently generated by computer 
“preprocessors”, and the even more voluminous output can be effectively 
analysed and presented by “postprocessors”. Hence a larger problem does not 
involve undue additional effort for users. 

(4) Problems of size and complexity hitherto unimaginable and infeasible can be 
handled by FEM, enabling analysts to extend their investigations into fresh 
areas, and inspiring designers to create new forms and new solutions. 

(5) Where formerly only a few alternatives could be examined, with FEM quite a 
large number of possible solutions could be tested, resulting in optimal solutions. 
However, certain disadvantages and limitations of FEM should also be 

recognised: 

(1) Every finite element is based on an assumed shape function expressing internal 
displacements as functions of nodal displacements. A certain element may give 
accurate answers for a particular type and location of support and loading, but 
inaccurate answers for another type and location. 

(2) Even with “well-behaved” elements, the solution is heavily dependent on the 
mesh, not only on the number of elements into which the region is divided, but 
also on their shape and arrangement. 
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(3) Modelling of the geometry, material properties, support conditions, and loading 
conditions, is very subjective and depends on analyst's judgement a lot. The 
same problem, when solved with the same computer program by different 
analysts, will often result in answers differing to a smaller or larger degree. 

(4) Precision of the output to a large number of significant digits 1s no guarantee of 
the accuracy of the solution. Even convergence with refinement of mesh is not 
absolute proof of the correctness of solution. 

(5) As finite elements become necessary only when other methods are not available 
or economically viable, in many cases there will be not many ways to confirm 
whether the solution is right or exact. 

Hence FEM needs considerable experience in the problem modelling and result 
interpretation aspects. In the wrong hands it can lead to wrong and even dangerous 
results. One must simply investigate the system with various mesh sizes and types, 
different modelling treatments, and under all possible material, support and loading 
conditions, and look for consistent and convergent behaviour. 

Before results from any analysis are used in real-life situations, they must be 
backed by justification of the model, evaluation of convergence, and an estimate of 
the accuracy, at least in relative terms. 

At the same time, there is no cause for undue alarm. Fortunately, in the last few 
decades, innumerable types of problems have been analysed by many researchers 
and application engineers around the world, and a massive amount of practical 
experience has been accumulated in applications of FEM. 

Most of the packages have built-in checks for catching obvious errors. Hence the 
chances of a user making a large error in applying FEM are very little. 

As long as users are aware of the factors that could influence the result and 
techniques for eliminating or at least minimising the errors, results from FEA should 
be reliable. 


1.8 Mathematical Formulation of Finite Element Method 


As already mentioned, there is no theory that can yield the force-displacement 
relationships for a general triangle or quadrilateral. We have to make some reasonable 
assumption on the behaviour of the element and then use the fact that as we 
increase the number of elements the answer converges, to lead us to a usable 
solution. 

In the stiffness approach, the assumption is usually made at the displacement 
level. This will ensure compatibility of displacements between elements but not 
continuity of stresses. Rarely and in special cases, stress distributions may be 
assumed and a flexibility analysis made. More often however, a hybrid approach is 
used, enforcing partial contmuity of both stress and deformation distributions across 
element boundaries. 

In this book, the pure displacement approach will be described, with comments 
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on hybrid elements where necessary. 

In the displacement or stiffness method, the internal displacements are assumed 
as functions of nodal displacements, and these functions are known as “Shape 
Functions”. It is on the basis of those shape functions that stiffness matrices and 
equivalent actions for element loads are developed. 


1.9 Shape Functions 


1.9, 1 General Requirements 


Let a point inside or on the boundary of a finite element displace by the amount 
{U} under some loading. The components of the internal displacement vector {U} 
will depend on the particular problem being analysed. 

For instance, in the case of plate bending it will be {0,, 6,, w}, the rotations 
about the x and y axes and the vertical (z) deflection. To illustrate, we shall use 2D 
plane analysis with deflection components i.e. x and y movements {u, v}, because it 
is the simplest and the most common case, and is easily extendable to the solid 
elements which are of the greatest importance to our topic. 

The general internal displacements {U} are related to nodal displacements {D}. 
Thus, for an element with nodes I, J,..., with nodal displacements (u;, v;), (u;, V;),°°°; 
the x and y displacements u and v at any point will be: 


u>= N,u;+ Nju;+ oer 
and v= Nivit Njvj+°"°. 


This relationship is called “Shape Function” and the relating matrix is designated 


[N]. 
{U}=[N]{D3 (1.13) 


Internal displacements of a point will depend on the location of the point, and 
hence shape functions will depend on the coordinates {x, y} of the point. 

The shape functions must be true for the entire element, in particular on its 
boundaries and at the nodes. Hence at the node I, where, x= x and y = y;, u must 
equal u;, and v must equal vj. This means that at I, N; must equal 1, and N; and N, 
must equal zero. 

Similarly, at J, Nj= 1, and N; and N,=0; and at K, N,=1, and N; and N;=0. 

Thus the shape function may be visualised as a surface which has ordinate of 1 
at the node to which it applies, and zero at the other nodes. Fig. 1.8 shows the 
shape function N; for the node I, for a rectangular element with (a) four nodes I, J, 
K, L, and, (b) eight nodes I, J, K, L, P, Q, R, S, illustrating the analogy. 

Obviously, the more nodes an element has, the more terms can be incorporated 
into the shape function, and the better the solution surface can approximate the real 
behaviour. 

For straight prismatic members and certain regular shapes, shape functions can 
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Fig. 1.8 Shape Functions for (a) 4-node Rectangle, (b) 8-node Rectangle 


be obtained from theoretical principles. But, for general shapes, some assumptions 
will have to be made, to be validated and refined subsequently. 


1.9. 2 Displacement Function Technique 


For certain simple cases such as a plane triangle, it is easier to express the 
displacements of a point in the element as functions of its position, rather than to 
assume a shape function directly. 

Then the shape functions may be developed through such “ Displacement 
Functions” as follows: 


Step 1: 
Write the displacements at any poimt as functions of its coordinates and 


undetermined coefficients, through an assumed “Displacement Function” [P], in the 
form: 


u=C,x+C, 
and v=CytC, or in general 
{U}=[P]{C} (1.14) 


with {C} being the vector of undetermined coefficients C,, C,,..., the number of the 
coefficients being equal to the number of DOF for the element, say m. 


Step 2: 
The displacement function Eq. (1.14) is true for the whole domain, and in particular 


at the nodes such as J, J, When we substitute the nodal coordinates and their 
displacements into Eq. (1.14), we get: 


u= Cx +C, 
and vi= CQ yitG, or in general 
tD}= [M]iC} (1.15) 
The [M] matrix is square, (n Xn) in size. 


Step 3: 
The coefficient vector {C} in Eq. (1.15) may be determined by inverting [M];: 
{C}=[M]' {D} (1.16) 
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Step 4: 


Upon substitution of Eq. (1.16) into Eq. (1.14), we get the internal displacements as 
functions of nodal displacements in the form: 


{U}= [P][M]' {D}= [N]{D} (1.17) 


defining the Shape Function as: [N]= [P][M]’ 

This is feasible where the [M] matrix can be inverted, and is of reasonable size 
for closed form inversion or economical computation. 

For all but the simplest of cases such as a plane triangle or rectangle, this 
approach is not practical and more general and direct methods have been discovered 
and developed. 


1.10 Element Stiffness Matrix 


Although we would use the triangular finite element 
IJK of Fig. 1. 9 to develop the general procedure for 
developing the finite element stiffness matrix, the 
technique will apply to other shapes with different 
DOF as well. 

As the entire element deforms, a typical interior 
point such as P = (x, y) moves to the position P’, by 
an amount {U}, with components {u, v}, as shown in Fig. 1. 9 Displaced Triangular 
Fig. 1.9. The rectangle around the point denotes an Finite Element 
infinitesimal element dxdy of the triangle. 

The general procedure for the determination of a stiffness matrix may be 
outlined as follows: 


1.10. 1 Shape Function 


Except for very simple cases, internal displacements {U} are written as an assumed 
function of node displacements {D}: 


{U}= [N]{D} (1.18) 


where, [N] is the shape function matrix. 

The individual shape function expressions will be functions of coordinates (x, 
Yi, Xm» Ym) Of the m nodes of the element and variable general point coordinates (x, 
y). That is, 


Ni= fi(xy), Nj= fy), (1.19) 
1.10. 2 Strain Influence Matrix 


Strains are derivative functions of displacements. For example, in the plane stress 
situation shown in Fig. 1.10, we have: 
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situation shown in Fig. 1.10, we have: 


au av av, au 
Bay EY Fye ME Pe ay ay 
Symbolically, the strain vector may be written as: 
{e}= [L][U] (1.20) 


. , ae d . 
where [L] consists only of the partial derivative operators such as 3x) and is called 


the “Operator Matrix”, implying that its terms operate on (rather than multiply) 
the terms of the displacement vector. 
From Eqs. (1.18) and (1.20), we get: 


te $= [LIN] {Ds= [B}iD} (1.21) 


where [B]= [L][N], representing the result of [L] operating on the shape function 
matrix [N]. The author will refer to the [B] matrix as the “Strain Influence Matrix” 
as it gives the strains when multiplied by the nodal displacements. 


1.10.3 Stress Influence Matrix 


Stresses are functions of strains by the generalised Hooke’s Law. In the plane stress 
situation of Fig. 1.10, we have: 


_ Elextvey) — _ Eve,tey) __Ety 
Ie ee ek Pe) 
(a) (b) ho (Cc) —- “Tr 
O. O- 
7 it 
dx 7 T 


dx dx 


Fig. 1.10. (a) Normal xstrain; (b) Normal y-strain; (c) Shear xy-strain 
Thus, the stress vector may be written, using Eq. (1.21): 
[o}= [E}te}= [EJ[B]{D}= [S]{D} (1.22) 


where [EF] is the elasticity matrix, involving the material properties E the Young's 
Modulus, v the Poisson's Ratio, and the product matrix [S] (= [E][B]) may be 
referred to as the “Stress Influence Matrix” because it gives the stresses when 
multiplied by the displacement vector. 


1.10. 4 External Virtual Work 


Now consider the element subjected to additional virtual (imaginary) displacements 
{D° } at the nodes, with the asterisk ( ) superscript indicating virtual quantities. 
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to be an elastically deformed region under actual loads {A}, and then an additional 
virtual displacement may be imposed on it, shifting it to the position indicated by 
the thin-line triangle J’K’. 

The external virtual work done by the actual actions {A} moving through the 
virtual displacements {D"} is given by: 


W. =D; A;+D, A.t+ = {D'}"{A} (1.23) 


the vector transpose {D"}' being necessary to convert the {D"} vector to a row 
matrix, thus enabling the actions and displacements to be multiplied and added, and 
hence to yield a scalar work quantity. 


1.10.5 Internal Virtual Work 


The nodal displacements will result in additional general internal virtual 
displacements {U°} and hence virtual strains {e'}. The actual stresses {o} in the 
element moving through the virtual strains {e"} will cause internal virtual work, that 
is to say, will store internal virtual strain energy. 

The same relationships that we developed for the real displacements and strains 
will also hold true for the virtual displacements also, and hence we may write: 


{U"}= [N]{D"} (1.13a) 
{e"}= [LI[N]{D }= [B]{D } (1.21a) 


Fig, 1.11 shows an element of unit sides and o 
unit thickness, so that stresses are equal in i 
magnitude to the forces, and strains are equal in | +» 
magnitude to the displacements. ee, 

For this unit cube, the virtual internal work X 
done (or strain energy stored) by the actual a ‘a 
stresses {o,, o,} moving through virtual strains 
{e. ,e, } is given by: y 


ponerse | 


Fig. 1.11 Real Stresses and 
Virtual Strains 


dWin= ex oxtey oy t= fe} {o} 


the transpose {e}" of the strain vector being 
necessary to convert it to a row matrix and yield 
the scalar work (or energy) quantity on multiplication. 

For an elemental volume of the body dV (= dxdydz), the virtual internal work is: 


dW? = fe"}" fo} dV (1.24) 


From Eq. (1.21a), we can write: {e"}'= {D"}"[B]' 
Substituting this into Eq. (1.24) and applying Eq. (1.22), for the entire finite 
element, the internal virtual work done, we get: 
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wi = | (DY BITE] [BYD}V 


VOL 


(D'}"¢| (BY EIBIav) (D} (1.25) 


VOL 


II 


1.10.6 Virtual Work Equation 


Total virtual work (W. —W;) must be zero. Hence from Eqs. (1.23) and (1.25) we 
get: 


(Dy tA} = (D'Y (| (BF LENBIav)(D} 


VOL 


or fA} = (| [BI TEIBIAV)(D} (1.26) 


VOL 


Eq. (1.26) may be symbolically written in the same form as Eq.(1.4): 
{A}=[K]{D} 


in which, [K], relating action and displacement vectors {A} and {D}, may be termed 
as the “Stiffness Matrix”, and may be determined for any element as: 


(K]= | (BY'[EIBIaV (1.27 
VOL 
Note that in the above, {A} and {D} are (m X1), and [K] is (m Xm), where,m is 
the number of DOF for the element, 6 in the case of the triangle considered. 


1.11 System Stiffness Matrix 


For an assemblage of elements, similar action-displacement relationships exist as for 
individual elements. Each action on the system is a linear function of the 
displacements at all the DOF connected to the DOF of the particular action. 

The basic approach to develop these relationships is the same, namely, to 
determine the actions at all the DOF in the system due to unit displacement at each 
DOF in turn. But this can become very tedious for large systems. Instead, a simpler 
method, based on equilibrium and compatibility at the nodes, is used. 

For illustration, consider the two elements P(= IJK) and Q(=IKL) of the 
machine part in Fig. 1.12(a). The DOF for all the elements of the system are 
sequentially numbered, and the eight DOF corresponding to the four nodes of the 
two elements P and Q are marked in Fig. 1.12(b). The curved surface of the pipe is 
now shown modelled as flat plate elements. 

IJKL is deformed to the shape IJ’K’L’ by actions(not shown) on this part as 
shown in Fig, 1.12(c), with displacement components {D,, D,,°**, Dg}. 

Fig. 1.13 shows elements P and Q with their separate nodal displacements and 
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Fig. 1. 12 (a) Machine Part; (b) Two Plane Elements IJK, IKL; (c) Displacements 


actions. The displacements at the common nodes I and K, namely D,, D,, Ds, and 
Ds, will be the same for both elements for compatibility, that is, to maintain 
continuity of the material. 

As already described, the stiffness matrices of the two elements can be 
determined and the nodal actions may be written in terms of the nodal 
displacements for each element. 

Thus, for instance, with superscripts within brackets identifying the elements: 


A= ky "D, i ki, D,+ k,;” D;+ kg Dyt kis“? Ds + kig Dg 
and, A, (Q) — ki, QD, ee k,@ D,+ k,,@ D,+ ke? D,.+ k 7? D,+ k, Ds 
Then, A,=A,“+A,@ 

ca (ky ies ki; S )D, 7 (+ ka )D,+ k,” D,+ k.? D, 
+ (kis ee k,;“ )Ds+ (kis? Ki” )Do+ Ky D+ Kis ~ Ds 
Write: A,= k,, Di + ky» D)+ ky3 D3+ ky Dat kis Ds + kig Det ki7 D7 + kig Dg 


A, 7 (Cc) 


(d) 


Fig. 1.13 (a) Displacements of IJK; (b) Displacements of IKL; (c) Actions for IJK;(d) 
Actions for IKL 


It is seen that the system stiffness coefficients at the common DOF 1], 2, 5, and 
6 are the (algebraic) sums of stiffness coefficients of the elements connected at these 
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DOFs. Thus, 
ky= ky ie ki wv ) kp= k,O+ | rae kis= ki+ Kis, kis= ke + kis 


Similarly writing the equilibrium equations for the forces along the other three 
common DOF 2, 5, and 6, we will get: 


a . = P oo P _— Pp 

k= ky,‘ + ky as ky = k»‘ ka kos= kos‘ ee Cee kos= kag Mt kee 
ens P — P _— P _ ly 

ks) = ks,‘ aks 5 ks. = ko + Ke | ks5= kgs‘ + ko ks5= kso\ + Keg 


kor = ke + ker, Keo= Ken + Ken, Kes= kes? + kes, Kes= Kes + Kes 


In fact, considering all the four elements connected at J 
the node K in Fig. 1.12(a), as shown in Fig. 1.14, the P 
stiffness coefficient ks, will not be just (ks + ks) as K 
given earlier, but will be as follows: ; / 

i O 


ks = ks + ks,“ + Ks; oa ks ° | : 
Fig. 1.14 Four Elements 


R and S being the elements to the left of elements P 
around K 


and Q. 

In general, a typical system stiffness coefficient in the system stiffness matrix 
[K,] may be written, summed up all the M elements connected along the DOFs I 
and J, 


ik, = Dk? (1.28) 
m= | 


More than the mathematical justification, it can be intuitively appreciated that 
for a node to displace along any of its DOF, it has to exert actions to overcome the 
stiffnesses of all the elements attached at that node. In other words, the stiffness at 
a node is the sum of the contributions from all the elements meeting at the node. 

The technique used in the computer assembly of the stiffness matrix for a 
system with n DOF from m individual member stiffness matrices may be now 
summarised as follows, with reference to Table 1.1 for the two elements P and Q: 
(1) Set up an (empty) system stiffness matrix [K,], (n Xn) in size; 

(2) Determine the stiffness matrices [K,,]| of all the m elements, identified with their 
appropriate system DOF; 

(3) Add each term of the stiffness matrix of every element into the row and column 
of the system stiffness matrix, according to the system row and column DOF. 

It schematically tabulates the two (6 6) stiffness matrices of the two elements 
P and Q in Table 1.1(Top), as well as the (8 <8) system stiffness matrix assembled 
from them in Table 1.1(Bottom). For ease of presentation, ki” and kK? have been 
shown as p,; and q;, respectively. 

Blank cells such as at the intersection of rows and columns for DOF (3,4) and (7, 
8) imply that the nodes corresponding to them, J and L in this case, are not directly 
connected. 

Thus, the system stiffness equation can be developed as: 
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{As }= [K.]{D, } (1.29) 


Table 1. 1 Element Stiffaess Matrices, (Bottom) System Stiffness Matrix 
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As the system stiffness matrix [K,] is simply the superposition of the 
component member stiffness matrices, all the characteristics of member or element 
stiffness matrices, listed in Section 1.5.2 in Chapter 1, carry over into the system 
stiffness matrix. In particular, the system stiffness matrix at this point is singular, 
and cannot be inverted. 
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1.12 Equivalent Actions Due to Element Loads 


The matrix formulation and solution can be only applied to express the effects of 
actions at the nodes. But in practical problems, the loads which cause deformations 
and consequent strains and stresses, may be on or inside the boundary of the 
elements. Such non-nodal loads are referred to as “Element Loads”. Examples are: 
(1) Concentrated action A,, applied on the boundary or within an element; 
(2) Distributed “Tractions” q, along the edge or surface of an element; 
(3) Distributed “Body Force” g, acting on all the particles of the element, as gravity, 
magnetic effect, or centrifugal action; and 
(4) Distributed “Initial Strain” ¢,, due to temperature change, force fit of slightly 
larger or smaller component parts, etc. 
The set of nodal actions that will cause the same deformations as each of these 
element loads, and hence may be substituted in their place, are referred to as 
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“Equivalent Loads”. 

We have already used the term in Chapter 1, during the discussion of the effects 
of uniformly distributed on fixed-ended beams in Section 1.5.3. 

We shall invoke the principle of zero 
virtual work to derive expressions for the 
equivalent loads. In Fig. 1.15, the IJK is 
defined as an elastically deformed shape of 
the finite element. Consider an additional 
virtual displacement {D"} imposed on each 
node. The corresponding equivalent nodal 
actions are {A,}. 

These will become are the equivalent nodal 
actions when the virtual work done by these 
actions equals that by the applied element 
loads. Note that each of the actions and 
displacements which have been identified with 
the node subscripts i, j, and k, will have the usual components along the coordinate 
directions. 


Fig. 1.15 Virtual Nodal Displacements 


1. 12. 1 Concentrated Action inside Element 


Consider the element IJK loaded with a concentrated 
force {C}, at the point P (x y), as in Fig, 1.16. Note 
that P may also be on the boundary. 

An additional virtual displacement described 
before is given to the element and the poimt P 
correspondingly moves to the position P’ by the ] 
amount {U"}= {u’ , v'}. 

External virtual work done by the applied action 
{C}, noting that {U"}=[Np]{Dz}, 


W.= {U"}"{C}= {D.}" [Ne ]’ {C} (1.30) 


Fig. 1. 16 Concentrated Force 
inside Element 


Then replace the applied force {C} by equivalent actions {A,}= {Aa, Aw,***} at 
the nodes (similar to {A.} of Fig. 1.15 so as to cause the same displacements at the 
nodes as {C} would. 

The virtual work done by the equivalent nodal loads through the virtual nodal 
displacements, 


W.= {D:}" {Ac} (1.31) 
Equating the two values of virtual work from Eq. (1.30) and Eq. (1.31), we get: 
{A.}= [Np]' {C} (1.32) 


Thus, the equivalent load vector {A,} may be computed with the shape function 
[N] evaluated at the location of {C}, and for the known value of {C}. 
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1.12.2 Traction on Edge of Element 


Let the element IJK be loaded with a traction 
{q}= {q., q,} per unit length, along the edge JJ, 
as shown in Fig.1. 17. Consider an element dS 
of the edge as marked with a thicker line. 

The element is given an additional virtual 
displacement {D'} at the nodes I, J, K, and 
then move to the position shown by the broken 
lines. The edge element dS correspondingly 
moves by the amount {U"}= {u' , v}. 

The applied traction {q} is replaced by 
equivalent actions {A,}= {Ag, Agw,***} at the 
nodes (similar to {A.} of Fig. 1.15) so as to 
cause the same displacements at the nodes as 
{g} would. 

External virtual work done by the applied traction {q}on the edge element dS, 
noting that {U°}=[N]{D'}, 


dW,= {U'}" {q}dS= {D'}"[N]' {q}dS 


Fig. 1.17 Uniform Traction 
along One Edge 


For the entire edge, external virtual work done by the traction is: 


w.= (Dy | INF (aids (1:33) 
[+] 
The integration to extend over the entire edge IJ subjected to the traction. 
As in the case of the concentrated action, external virtual work done by the 
equivalent actions {A,}, 


W,= {D }" {Aj} (1.34) 


Equating W, and W,, the two values of virtual work from Eqs. (1.33) and (1.34), 
we get: 


A.= | INI {ads (1:35) 


EDGE 


1.12.3 Body Force over the Element 


The element IJK is loaded with a body force {g} per unit volume, {g,, g,} in this 
plane case, over the entire element, as shown in Fig. 1.18. 

Consider an infinitesimal element dx dy of the finite element as marked. 
Elemental volume dV = tdxdy, t being the element thickness in this plane case. 

The element is given an additional virtual displacement {D"} at the nodes I, J, K, 
and then move to the position shown by the broken lines. The area element 
correspondingly moves by the amount {U"}= {u' , v}. 
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The applied body force {g} is replaced by 
equivalent actions {A,}= {A,, Ag} at the nodes 
(similar to {A,} of Fig. 1.15) so as to cause the 
same displacements at the nodes as {g} would. 

External virtual work done by the applied 
traction {g}on the volume element dV, noting 
that {U"}=[N]{D°}, 

dW.= {U'}" {g} dV 
= {D'}"[N]' {g} dv 

For the entire element, external virtual work 

done by the body force 1s: 


Fig. 1. 18 Body Force on Element 


W.= (Dy | (NI (gav (1:36) 
VOL 
the integration to extend over the entire area (and hence volume) of the element. 
As before, external virtual work done by the equivalent actions {A,}, 


W,= {D'}" {A;} (1.37) 


Equating W, and W,, the two values of virtual work from Eqs. (1.36) and (1.37), 
we get: 


A, = | [N]’ {g}dV (1.38) 


VOL 


1. 12. 4 Initial Strains in the Element 


Initial strains may be introduced to a body 
due to temperature change, the necessity to 
force fit a component into place, and so on. 

For example, when the temperature of a 
body changes by T, it experiences normal 
strains of a.T , where a is the thermal 
expansion coefficient of the material. No 
shear strains are developed. 

Then the element IJK is loaded with 
initial strain {eg}= {€, €y0, Yyo}, over the 
entire element, as shown in Fig, 1.19. 

Consider an infinitesimal element dx.dy of the finite element as marked. As 
before, elemental volume for the plane case, dV = tdxdy. 

The element is now given an additional virtual displacement {D"} at the nodes I 
J, K, to the position shown by the broken lines. 

The applied initial strain {e)} is replaced by equivalent actions {A,}= {A.,, Ax} 
at the nodes (similar to {A.} of Fig. 1.15) so as to cause the same displacements at 


Fig. 1.19 Initial Strains in Element 


> 
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the nodes as {e)} would. 

In the case of initial strain, it must be remembered that only the elastic 
component of the strain (total strain minus the initial strain) would cause stresses in 
the material. 

Hence Eq.(1.22) must be modified as: 


[o}= [E] fe ~€03= LEN[B]{D}—[E] {eo} (1.39) 
Then, the expression for the internal work done, Eq. (1.25) gets modified as: 


Wi = {D" JKC | (BITEI (BI avy} 


YOL 


- | [BI TElfe0} av) (1.40) 


VOL 


The virtual work equation leading to the stiffness relation Eq. (1.26) is modified 
to: 


(A}=( | [BITE] (B] av(D}- | (BI TEI{e0} av 


VOL VOL 


=[K]{D} ~ | [BI'[E}fe0} dv (141 


The second term on the right-hand side represents the effect of initial strain and 
has the dimensions of force, or action in general. It may be designated as the 
equivalent load {A,} corresponding to initial strain, and shifted to the left-hand side, 
to give: 

{An }+ {A }= [K]{D} (1.42) 
where {A,} represents the actions directly applied at the nodes, and 


(AJ= | [BI LEIeo} av (1.43) 


VOL 
Incidentally, [E]{e.} may be a directly applied “Initial Stress Vector” {o)} also. 
1.12.5 Total Action Vector 


Thus for an element subjected to all the four element loads described in the 
preceding sections, the total action vector will be: 


{A}= {An}+ {Apt (Ag}t tAg3t {Ac} (1.44) 


where {A,} equels concentrated actions applied at nodes, and {A.}, (A,}, {A,}, and 
{A. }are known from Eqs. (1.32), (1.35), (1.38) and (1.44). 

Eq. (1.44) is evaluated for all the elements, and then summed up to give the 
system action vector {A,}. The stiffness equation will continue to be: 
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tAs $= [Ks ]{D; 5 


with stiffness matrix [K,] being the superposition of individual element stiffness 
matrices given by Eq. (1.27). 

Further solution will proceed according to the methodology described in Section 
1.6 of Chapter 1, to yield displacements and support reactions. 


1. 13 Stresses and Strains 


Apart from the deformations of the body analysed, the critical information always 
needed is the ability of the body to resist the applied loads without yielding or 
breaking. 

In matrix analysis of structures, the end results of matrix analysis are the 
internal axial and shear forces and bending moments in the members. In finite 
element analysis of continua, the end results are strains and stresses. 

Strains are simply obtained from Eq. (1.21), with the [B] matrices recalled from 
storage, if stored earlier. If we are only interested in a few critical locations, it may 
be cheaper to re-compute the [B] matrices for the specific elements. 

Stresses can likewise be obtained from the strains using Eq.(1.22). Storage or re- 
computation of the stress influence matrix [S] is not necessary. 

The numerical output from finite element analysis will be overwhelming, too 
much to understand and digest. These days, with so many automated graphical 
facilities available on the computer, most packages have modules to present the 
results in graphical form. Deformed shapes, with displacements increased hundred- 
to thousand-fold to make them visibly meaningful, and strain and stress 
distributions displayed as contours with lines and colour bands depicting various 
value ranges, are excellent tools for understanding. 


1. 14 Stiffness Matrices for Various Elements 


In many standard texts on the finite element method, at this point, actual stiffness 
matrices are usually developed and/or presented for various shapes and types of 
elements by application of the principles described in the preceding sections. 

This would include: Plane stress, plane strain, axisymmetric, solid, plate 
bending, and shell elements, with various DOFs; and refinements such as 
isoparametric formulation, varying orders of integration, etc. 

However, in this book, such development or even presentation of stiffness 
matrices will not serve such purposes. If readers wish to learn more about these 
specific topics, they may consult standard references such as Przemieniecki, J.S. and 
Zienkiewicz, O.C. 

In each of these, in spite of the mathematical rigour of the actual development of 
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the stiffness matrix, the starting point is the shape function, which is an assumption 
and this is why the finite element method requires caution in its practical 
applications. 


1.15 Critical Factors in Finite Element Computer Analysis 


Computers are used for to sdve problems with a large amount of calculation; finite 
element method is used to solve complex problems. Their combination leads to some 
special factors that must be considered and resolved, such as the following; 


Computer modelling 


The process of transforming and transferring details of an actual system into a form 
that the computer can recognise and work upon is called “computer modelling” or 
“idealisation”. It must be remembered that the computer will never know or handle 
the object itself. It can only manipulate numbers representing the various parts of 
the object and the actions which it is subjected to. 

As such, the way users translate the object into mput information will be the 
most critical aspect of the analysis. In certain cases, the modelling may be obvious 
and suggested by the configuration of the object itself. But in objects of continuous 
solids of a general irregular shape, much care and planning become important. 
Validation of the model may be critical to the acceptability of the results. 

Input preparation | 

Large problems require large amounts of data, in terms of the coordinates of large 
numbers of nodes and the incidences (i.e. node connectivities) of large numbers of 
elements, as well as other items relating to material properties, boundary conditions, 
loading conditions, etc. Special means must be taken to produce and transfer them to 
the computer quickly and without error. Programs used for this purpose are called 
“Preprocessors”’. 

Out put evaluation 

Similarly, the results from the analysis will also be massive, which must be 
evaluated carefully to identify the magnitude, nature and location of worst effects, 
again quickly and correctly. In finite element analysis of unusual problems, ways 
must be provided for users to view summaries of critical results. Where the object is 
complex and assumptions (as in support conditions) have been made, some means of 
bounding the magnitude of possible error must be provided. Programs used for this 
purpose are called “Postprocessors”. 

As a case study on all three aspects cited above, the author cites his paper on 
the finite element analysis for the stress concentration in a steel pressure vessel. 


1. 16 Modelling Considerations 


Most finite element problems are unique and deserve special attention in modelling, 
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Even experienced analysts who investigate a finite element problem using the same 
software, may come up with varying results sometimes differing by orders of 
magnitude! 

Except for the simplest and most conventional cases, modelling of a structure or 
system for computer analysis is a skill involving subjective decisions and 
interpretations. The main fact to be remembered is that the computer deals only 
with and in numbers, without the implication of any physical quantity or real 
entity. 

Hence two important steps must be included in the modelling/analysis cycle, 
namely : 

(1) Model validation, and 
(2) Result verification. 

With computers, engineers and scientists can analyse larger and more 
complicated systems, and include more of the surrounding domain, more of the 
supporting media, and so on, than before the advent of computers. 

Under these circumstances, idealisation becomes even more important than when 
the structures were simpler. Furthermore, intelligent assumptions and decisions will 
still have to be made on the following questions: 

(1) Should a 2D analysis or 3D analysis be done? 
(2) Should every little detail be included, and if not, what can be omitted? 
(3) How much of the supporting medium should be included? 

The question of 2D vs. 3D is a subjective one. With today’s powerful 
computers, users tend to go to 3D as a rule. However, while it may seem more 
realistic, there are two counter-arguments: 

(1) Is 3D really necessary all the time? 
(2) Does it involve extra information, further assumptions, etc. in the input? 

The answer is often “No” to the first, and “Yes” to the second one. There is 
considerable published material on this topic to convince one to use discretion and 
experience before launching a full-blown 3D effort. 

The author has occasionally run various cases of a particular problem by both 
2D and 3D, and extracted a “Correlation Factor” between 2D and 3D (for specific 
results), and later run only 2D cases of the same problem, applying the correlation 
factor to predict the more realistic 3D result. 

One difficulty which is unique to finite element analysis is that the magnitude of 
the “discretisation error” arising from the division of the domain into elements is 
often unpredictable, and needs independent validation. Errors may also occure due 
to wrong choice of type of element, wrong modelling of the boundary or loading 
conditions, etc. 

The best protection against making a wrong choice of mesh or input parameters 
is experience developed by long use, or knowledge borrowed from experts, through 
publications or otherwise. 

Lacking these, users may need to adopt one or more of the following methods to 
evaluate the validity of the computed findings and extract dependable results from 
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finite element analysis: 

(1) Physical testing some typical continua of the type and general dimensions 
analysed and comparing the experimental and computed results. But test 
specimens and tests are often too costly, involve delays, and may or may not 
conform to the object and conditions of the actual problem. 

Furthermore, most of the time, finite element analyses are resorted to only 

when the other classical or matrix methods cannot be applied, and as such, it 

may not be feasible to apply experimental methods to them. 

(2) Making the finite element model self-validating as follows: 

1) Sensitivity studies must be made on computer models including various 
parameters and their combinations likely to influence the object's behaviour, 
and their results evaluated to ascertain which of the parameters will be more 
critical and which less. 

If this is not done, important parameters may be left out by oversight or 
assumed insignificance. Alternatively, any and every parameters may be 
included, over-burdening the computer input, processing, and evaluation. 

2) Convergence studies must be conducted on at least two meshes to obtain an 
estimate of the discrepancy in the critical stresses. 

If the difference in results from two different meshes is beyond tolerable 
percentage, a third (perhaps even a fourth) finer mesh analysis will be needed 
to check the convergence of the results to some asymptotic value. In the 
author's experience, most problems give rise to wide variations for the same 
problem from different meshes. 

In this case, the results for critical stresses or displacements from three or 
more meshes must be utilised to extrapolate an improved (predicted) result 
for the continuum, as described later, and as demonstrated in the author’s 
paper (Ref.3.1), and summarised in a later section. 

(3) In many cases, solution to a similar (simpler) configuration of the subject case 
may be available by other methods, or from other sources. Such problems may 
be analysed by finite elements, and the type and density of mesh that best 
predicts the known values may be adopted for the complex case on hand. 

(4) Occasionally, the configuration of the object may be close to a problem which is 
amenable to classical or approximate solutions. In such cases, an analysis may 
provide a rough check on the validity of the finite element solution. However, 
with complex structural forms and combinations commonly analysed by finite 
element method, approximate solutions may be far away from the exact solution, 
and comparison with computer solutions in such cases may be misleading, 

Then, “self-validation” by running a few more cases of the same problem 
with slight modifications of supports, connectivities, element idealisations, etc. 
may be quite informative. 

Once a model has been validated, parameter studies wherein significant 

parameters governing the problem are varied between their maximum and minimum 

practical values (and preferably even a middle value), and the results for significant 
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variables evaluated by regression analysis, to obtain good prediction equations for 
the continua. 


1.17 Asce Guidelines 


A publication of the American Society of Civil Engineers gives some guidelines for 

finite element idealisation. Some relevant ones are listed hereunder: 

(1) Finite element gridworks should be uniform where this is practical. However, 
nonuniformity is often necessary to achieve refinement of the gridwork in 
regions of rapid changes in loading and/or behaviour. 

(2) Quadrilateral finite elements should be used except where triangular shapes are 
necessary to effect grid refinement or accommodate irregularity mm geometry and/ 
or loading. 

(3) Analyses to obtain displacement (or force) behaviour do not require as refined a 
gridwork as is required to obtain strain (or stress) behaviour. 

(4) Analyses including nonlinearity in geometry or material usually require a more 
refined idealization than linear analyses. 

(5) Where possible, gridpomts should be numbered from regions of large displacement 
toward regions of small displacement (e.g. from tip toward root for a cantilever 
beam). 

(6) Monotonic convergence characteristics are seldom guaranteed with respect to 
behaviour at a point; an exception may be the displacement under a single 
concentrated load. 

(7) When flat elements are used in the polygonal approximation of curved shell 
surfaces, rotation about the normal to the curved surface should usually be 
retained as a degree of freedom. However, the retention of such degrees of 
freedom may yield an illconditioned stiffness matrix. If these degrees of 
freedom are suppressed, it is important to check that the resulting reactions are 
acceptably small. 

(8) The use of constraint equations as a means of imposing special connections and 
boundary conditions may increase the bandwidth substantially. 

(9) The use of constant strain triangle finite elements should be avoided wherever 
possible. 

(10) The recommended first checks on predicted behaviour are equilibrium checks 

using loads, reactions and forces. 

(11) Aspect Ratio (length/height) limits vary from element to element; however, a 
general guideline is to keep element aspect Ratios below 10 for deformation 
analysis and below 5 for stress analysis. 

(12) Curvature of sides of iso-parametric elements should be employed only to 
represent curved boundary geometry or curvature of stress iso-clinics (where 
this can be reasonably inferred from the geometry and loading). 

(13) Non-uniform distribution of intermediate nodes on the sides of higher order 
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elements should be avoided wherever possible. 

(14) With higher order elements, special care must be exercised in defining boundary 
conditions because the solutions for some elements are sensitive to this choice. 

(15) Elements of extreme shape (high aspect Ratio rectangles, large angled triangles, 
re-entrant corner pentagons) are unfavorable with respect to manipulation error 
and should be avoided. 

(16) If the structure and loading possess symmetry, advantage should be taken into 
account of this fact. 

(17) For anisotropic materials, the subscript on Poisson's Ratio must be clearly 
defined. In addition, the theoretical limits on v, E and C should be checked. 
(18) For isotropic materials in “compressible solids” programs, the value of Poisson’s 

Ratio must not be too close to 0.5. 

(19) Rapid changes in element size should be minimized. 

(20) Do not try to fake internal openings by using a very soft (i.e. very low E) 
material unless all the “free” nodes in the opening are restrained. (This is used 
in conjunction with some mesh generators which do not permit internal 
openings.) 

(21) When modelling reduced portion of the structure for evaluation of local effects, 
make the boundary of this local area far enough away from the area of interest 
to avoid boundary effects on the information at the area of interest. 

(22) Connecting high order elements to low order elements can cause stress 
irregularities which should either be smoothed or ignored. 

(23) For complex structures consider a “staged” analysis. Analyze the entire 
structure with a coarse mesh and then study the areas of interest with a fine 
mesh, using as boundary loads the results of the coarse analysis. 

(24) If flat shell elements are used to model a curved shell, the pressure load should 
be applied as concentrated loads at the nodal points. 

(25) In modelling long cylindrical shells, the model should include at least three 
wavelengths on either side of any perturbation in geometry and/or loading, 
(26) Two elements through the wall thickness of a thick-walled pressure vessel will 
usually give reasonable results with incompatible displacement elements. 

(27) Solutions for structures with uniform stiffness distribution are more accurate 
than those for structures with stiff and flexible regions. 

(28) A boundary layer of width equal to plate thickness usually exists along the 
edges of shells. There are high stress gradients in the boundary layers. To 
properly reproduce these stress gradients, there should be a double row of 
elements along the edge, the overall width of the row being equal to the 
boundary layer width. 

(29) A boundary layer of width t may exist along the edge of a plate which 
undergoes shear deformation, in addition to bending deformation. Remarks 
similar to Guideline 28 apply. 

(30) If the solution routine is band-width dependent, the band width should be 
minimized. Likewise, if a solution routine is dependent on the element 
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numbering scheme, an optimum sequence should be used. 
(31) For convergence studies, one must refine (i.e. sub-divide) the original mesh, 
otherwise a new approximating sequence is begun. 

Most of the guidelines must be taken only as suggestions or recommendations, 
but not as rules, because with finite elements, decisions must be made based on 
merits of each case, unless it happens to be almost identical with others already 
successfully analysed earlier. 


1.18 Preprocessors and Postprocessors 


1. 18. 1 Preprocessors 


Fig. 1.20 depicts a machine part, a continuous solid that has been divided into a 
mesh of 188 solid elements, connected at 396 nodes. This is given only for 
illustrative purposes. Actual finite element models would be much more complex 
and have thousands of nodes and elements. 


Ete ita 


Fig. 1.20 Machine Part Discretised into Finite Elements 


As listed in Section 1.5, the input normally includes: 1) Geometry, 2) topology, 
3) material properties, 4) element types and cross-sectional properties, 5) boundary 
conditions, and 6) loading conditions. 

Each item would typically involve an average of five to ten items of data for 
each node and element. The total could run to millions of numbers, every one of 
which had to be right for the analysis program to work, and then for the answers to 
be correct. 

Even if a large number of people were employed to generate these numbers, it 
would be a monumental if not impossible task to detect and eliminate all the errors 
that would inevitably creep in. 
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So computer programs for special purpose are used to generate and display the 
data for matrix structural and finite element programs for the analysis of different 
objects. These programs are called “Preprocessors”. 

In due course, computer graphics offered powerful display and manipulative 
capabilities. Soon after computer analysis package developers added preprocessor 
modules to display and edit the input and generated data on the monitor. 

Now the analyst could literally see his model and confirm the correctness of his 
data, or locate and correct the errors, and then proceed with the analysis without 
delay of time or waste of computation. 

The automatic generation of data in regions of regular geometry and topology, 
selective data display, zooming and rotation facilities, text annotation and colour 
coding, and on-screen editing capabilities, have now freed the analyst from struggling 
with the data preparation, checking, and correction. 

One might almost say that without these preprocessors, practical application of 
the matrix and finite element analysis would not have developed and been applied to 
the extent it has been applied and it is being applied today. 

Current preprocessors may be categorised into the following broad groups: 

(1) Numerical routines to interpolate omitted node coordinates and element incidences at 
equal or specified (proportionate) intervals. 

(2) Graphic input modules by which users can “create” and sub-divide his structure. 

(3) Menu-driven text help (input instructions) mode to guide users into developing 
the input commands. 

(4) “ Wizards”, namely user-friendly, icon-based, menu-driven guidelines for 
automatic generation of data required for the analysis, based on a few input 
parameters from the user, mostly by mouse-click selection. Many assumptions 
are built in to minimise input data. 


1.18.2 Postprocessors 


If the situation of the large input was bad, the problem with the voluminous output 
became worse. For each node and each element, large number of displacements, 
forces and reactions, moments and torques, stresses and strains, etc. were computed. 

Printed output for a single analysis could fill boxes, and it was impossible for 
analysts to go over all the numbers and keep track of the changes, to find the critical 
values necessary for design. It was too easy to get lost in the maze of numbers and 
miss the big picture. 

So analysts sought selective output in critical regions. Large number of 
draftsmen toiled to plot critical quantities across specified sections and lines. Again, 
computer programs were written to digest the massive output into manageable and 
understandable summaries. These came to be called “Postprocessors”. 

Here too, computer graphics came to the rescue, to even greater effect than with 
the input. Results could be superposed over the image of the object, also selectively, 
zoomed up in specific regions, and displayed in desired orientations. Images of the 
distorted object, force and moment diagrams for articulated structures, and stress 
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and strain distributions for continuous solids could all be viewed on the click of a 
button. 

This feature enabled immediate validation of the computer model and interactive 
analysis for changed configurations and for different boundary and loading 
conditions, eventually leading to optimum design. 


Fig. 1.21 Spanner (a) Deformed Shape, and (b) Stress Contours 


Fig. 1.21 depicts the computer “Scope” display from GTSTRUDL (Georgia 
Tech STRUctural Design Language) an integrated structural and finite element 
analysis package, for the analysis of a spanner tightening a hexagonal nut. 

The left sketch shows the original and deformed shapes. The right sketch 
displays the stress distribution contours. 

This is just an illustrative example, as real-life situations would be much bigger 
in data size and complexity. 


1,19 Support Modelling 


In earlier times, any component being analysed was assumed to rest on a firm 
foundation, often fixed to the base. In reality, however, there is a large interaction 
between the components and their support, the magnitude of which may be quite 
significant. Modern computerised methods enable the designer to take account of 
these effects. 

It is very important to decide whether a node is constrained in all its DOFs 
(namely “ fixed”), or constrained on certain DOF (“pinned”, “sliding”, etc.), or not 
constrained at all(“free”). In axi-symmetric solids or shells and in symmetric halves 
of components, the support conditions may not be obvious, and wrong choices will 
lead to very erroneous results. 

The support medium may also be modelled by finite elements with actual 
properties being input. The extent of the support medium should be at least 2 to 3 
times the size of the component being analysed, and would involve quite a bit of 
discretisation. Alternatively and with sufficient accuracy, they are usually modelled 
as springs under the component. 


1.20 Improvement of Results 


The result from analysis with a single mesh, unless substantiated by other evidence, 
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is really good only for general qualitative interpretation, and may not be sufficiently 

accurate for dependable quantitative conclusions. 

It is well known as one of the strengths of the finite element method, that for 
well defined problems, the results converge asymptotically as the mesh is refined, 
that is as n, the number of DOF (or the number of nodes or elements) is increased. 

Many studies have tried to develop a means to extrapolate the results from two 
or more meshes to the asymptotic result, assumed to be “exact”, corresponding to 
an infinite number of elements corresponding to the actual continuum. 

One conclusion is that when a mesh is progressively divided so that a mesh is 
developed from another mesh by sub-dividing the coarser equally in all coordinate 
directions, then the error (for a certain class of elements and for square shapes) is 
the order of h’, where h is the element size parameter. Note that the element size is 
(roughly) proportional to (1/n). 

Thus, if a rectangle 6 cm by 4 cm is first analysed with square elements of 2 cm 
side, and then analysed with square elements of | cm side (obtained by sub-dividing 
the original mesh equally in both directions), then the error is reduced to (1 cm/2 
cm) or one-fourth of the starting value. 

This observation has been utilised to extrapolate known results for two or more 
meshes to the continuum by what is known as “Richardson’s Extrapolation”, which 
has been widely used in finite difference applications. 

However, in the general case, square (or cube shaped) elements are not the norm, 
where mesh is often not uniform but graded, and it is very difficult to get mesh 
sequences increasing 4 times for 2D or 8 times for 3D for every successive 
refinement. 

Hence in general situations not conforming to the theoretical requirements of 
mesh variations, a strategy successfully adopted by the author is a modification of 
Richardson’s Extrapolation as follows: 

(1) Solve the problem by three or more meshes with DOF numbering n,, m, m,..., 
not necessarily in any sequence of sub-division. (Note that the technique must 
not be applied to results from two meshes!) 

(2) Confirm that corresponding values of any critical result, say s:, S:, 53,...tend to 
converge. Fig. 1.22(a) shows the plot of raw data. Note that at this stage, 
attempting to estimate the asymptotic value may be very subjective, as indicated 
by the two possible curves shown by full and broken lines in the nets both 
passing through the three points. 

(3) Plot s;, S2, S3,...against (1/n,)’, (I/m)’, (1/y).... 

(4) If the three points happen to be in a straight line (or at least nearly so) as shown 
by the full line in Fig. 1.22(b), then extend the line through them backwards to 
the vertical axis, to obtain the best estimate for continuum from the three 
meshes. 

(5) In general, the three points will not lie on a straight line, but in a curve as shown 
by the broken line in Fig. 1.22(b). By trial and error, find the value of the 
exponent m in (1/n)" for which the three points lie in a (nearly) straight line. 
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Result (a) Result (b) 
Asymptote *™~ Best Estimate 


No. of DOF 


Wy nA, 


Fig. 1,22 Results from Three Meshes for (a) Linear Plot, and (b) Exponential Plot 


Then extend the line backwards to obtain the best estimate. 
This procedure is simply a practical way to get a best estimate from three (or 


more) sets of results on hand, and not a theoretically valid method to obtain the 
asymptotic value. Furthermore, such extrapolation and the estimated value are valid 
only locally and for the specific critical value plotted. The more the number of 
meshes, the better the estimate of the asymptotic value. 
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2.1 History of Dental Implants 


Dentistry aims at replacing missing teeth since it was first recognized as a 
profession. For centuries, dental practitioners have relied on their own skills and 
various artifacts to develop esthetic and functional alternatives to minimize sequelae 
that occur as a result of edentulism. Partial, complete, fixed, or removable dentures 
are by far the most commonly used forms of tooth replacement applied. In other 
words, these devices have been incorporated into the oral cavity anchored on either 
remaining teeth and/or other anatomical structures. Only scarce archeological reports 
have demonstrated attempts of mcorporating prosthetic devices into the jaws as 
more natural and functional replacements. However, predictability of these methods 
was not achieved until recently. (Table 2.1) 


Table 2, 1 Advantages & Disadvantages of Implant Therapy 


Advantages Disadvantages 
Highly predictable Surgical intervention 
Preservation of osseous structures Higher costs 
Natural masticatory function Esthetics 


In 1952, in a modestly appointed laboratory in the university town of Lund, 
Sweden, Professor Per-Ingvar Branemark had an extremely fortunate accident, what 
most scientists call serendipity. Much to his irritation, Dr. Branemark observed that 
it was virtually impossible to recover any of the bone-anchored titanium microscopes he 
used in his orthopedic research. The titantum had apparently bonded irreversibly to 
living bone tissue, an observation which contradicted contemporary scientific 
theory. Metal implants in contact with living bone were thought to be embedded in 
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soft tissue, a term known as fibro-osseous integration. This form of metal-to~bone 
integration was never capable of achieving strong anchorage to bone. Dr. Branemark 
subsequently duplicated his findings and demonstrated that under carefully 
controlled conditions, titanium could be structurally integrated into living bone.’ 
This titaniunrto-bone integration was found to be extremely rigid, in contrast to the 
previously known concept of fibro-osseous integration. Dr. Branemark named this 
phenomenon. Osseointegration has been described as a direct structural and functional 
connection between living bone and the surface of a load-bearing implant. This 
concept was later modified with more advanced tmaging techniques. Scanning electro- 
microscopy revealed that a thin proteoglycan layer exists between titanium and bone. 

Soon after presenting osseointegration, Dr. Branemark realized that his discovery 
would have uncountable applications beyond the field of orthopedics. The most 
significant application would probably be anchorage to intra-oral prosthesis. 
Leonard Linkow was the first to utilize titanium in intra-oral rehabilitation. He 
reported utilization of titanium blade-form intra-osseous implants to provide 
stability and enhanced function to partial and complete dentures.’ These findings 
triggered developments of a series of studies to better illustrate the potential of 
osseointegration in treating edentulism.”* Later, Dr. Brdnemark revisited this 
concept by introducing hollow titanium screws for the same purpose.® An increased 
success rate, clinical applicability, and reduced rate of complications were obtained 
with this type of fixtures when compared to previously used blade-form implants.’ 
This screw has a very similar morphology and configuration with currently used 
endosseous implants. The hallmark paper of modern implant dentistry was 
published by Dr. Branemark in 1977 in which the first long-term longitudinal study 
in this field was reported. This 10-year report provided scientific data for the 
development of the currently used implantation protocols and implant systems.® 
The replacement of missing teeth with implant-supported prostheses has been a 
widely accepted treatment modality. Implants as a treatment option will be extensively 
used due to their high success rate and wide array of restorative options.” 


2.2 Phenomenon of Osseointegration 


The interactions between synthetic materials and host tissues require achievement of 
an optimal biological interface between the implant and the surrounding tissue. 
Adhesion molecule/cytokine interactions, specially osteopontin, an adhesion 
molecule as well as a cytokine, seem to play an important role during the wound 
healing process of dental implants’. An acute inflammatory response occurs after 
implant osteotomy preparation that includes the releasing and activation of a variety 
of cytokines and growth factors, which mediate initial events. Most of these factors 
come from platelets, and the others are released by PMNs and macrophages. As a 
result, mesenchymal cells migrate onto clots and synthesize a collagen network that 
becomes a scaffold for wound repairing. Local factor production determines the 
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quality of bone formation or the formation of fibrosis’. (Fig. 2.1) 
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Fig, 2.1 Sequence of Events During Osseointegration of Dental Implants 


The earliest events following implant placement are adsorption of water, ions 
and small biomolecules and are followed by adsorption and exchange of larger 
biomolecules such as proteins. The properties of the implant surface and those of 
individual proteins dictate the organization of the adsorbed protein layer, which in 
turn determines the cellular response to implant surface. 

Following adsorption of water, ions and small biomolecules, cell proliferation 
begins. Cell proliferation is modulated and influenced by TGF-81, PGE,, BMP-2, 
and 1,25-(OH), D; (vitamin D). The proliferation of undifferentiated mesenchymal 
cells is followed by differentiation of osteoblast-like cells. Appropriate mineral 
deposition is critical to ensure that the subsequent remodeling cycle on the implant 
surface will be optimal. Matrix production is necessary for the correct 
interrelationship between the apatite crystals and collagen’. Integrins are 
transmembrane proteins that mediate attachment to extracellular matrix proteins. 
Integrins as ECM receptors are believed to play a role in signal transduction and 
gene regulation, potentially influencing cellular migration, differentiation, and 
maturation. Cells are not actually attached to titanium itself, but to the organic 
material that adheres to the surface of that material. In this way, they form focal 
attachments through their integrin receptors, with arginine-glycine-aspartatic acid 
(RGD) sequences in proteins such as fibronectin and vitronectin”. 

The formation of a TiO, layer is considered to be the key to ideal interaction 
between titanium and bone’’. The result of dissociative adsorption of water on TiO, 
is a hydroxylated surface with two different types of hydroxyl groups. TiQ, is an 
amphoteric surface and has the ability to bind with both acids and bases". For cell 
attachment to be possible, the active sites of proteins have to be accessible to the 
corresponding sites on cell membrane, which means that globular proteins such as 
fibronectin have to be adsorbed in a non-denaturated state similar to their 
conformation in solution. The high biocompatibility of titanium can be explained 
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since the amorphous TiO, layer is a favorable surface for protein adsorption”. Four 
interrelated properties of biomaterials (chemical composition, surface energy, surface 
roughness, and surface topography) affect mesenchymal cells in vitro. Attachment, 
proliferation, metabolism, matrix synthesis, and differentiation of osteoblast-like cell 
lines and primary chondrocytes are sensitive to one or more of these properties”. 
Osteogenesis around implants is affected by the physical and chemical characteristics of the 
biomaterials used. The osteoprogenitor cells must migrate to the implant site and 
synthesize and secrete a mineralizable extracellular matrix. Because this is new bone 
formation, the mechanism by which the cells calcify their matrix involves extracellular 
organelles called matrix vesicles in a process termed “primary mineralization”*'. However, 
this process has shown to be delayed when compared to normal primary bone 
formation during bone healing. The delay in mineralization is compensated by an 
enhancement of primary mineralization by titanium”. 

The healing response of bone to the injury of implant insertion is similar to its 
response to fracture. If there is absolute stability between the bone and the implant, 
the injury to the bone will result in bone formation conforming to the surface 
configuration of the implant. There is first appositional bone growth and 
subsequently bone ingrowth into the interstices of the surface. The newly formed 
bone is not an uninterrupted layer of bone. If the implant is relatively stable after 
insertion, the healing response results in the implant becoming enveloped by densely 
organized fibrous tissue with occasional area of fibrocartilage, without any area of 
resorption in the surrounding bone. If the implant is unstable after insertion the 
healing response causes the implant to become surrounded by loosely organized 
fibrous tissue with active bone resorption occurring in the surrounding bone”. 

The process of ossification around implants is composed of three stages. At first 
incorporation by woven bone formation occurs, followed by bone remodeling and 
maturation. Woven bone is the primitive type of bone characterized by random 
orientation of collagen fibrils, numerous osteocytes, and low mineral density. Woven 
bone formation dominates the first 4-6 weeks after surgery”. Lamellar bone is the 
most elaborate bone tissue. Packing collagen fibrils into parallel layers with 
alternating course gives it the highest ultimate strength. Lamellar bone replaces 
woven bone by apposition”. Intense remodeling in lamellar bone adjacent to 
endosseous implants characterizes the last stage of the healing process. Intense 
remodeling begins at around 3 months and decreases after several weeks, but 
continues for the rest of life.” Remodeling starts with osteoclastic resorption, 
followed by lamellar bone deposition.” 

Little information is available on the bone-resorptive activity that occurs around 
implants during osseointegration. However, the role of osteoclasts, macrophages, and 
stromal cells in those catabolic processes associated with bone remodeling has been 
demonstrated.” The integrity of connective tissues surrounding dental implants may 
be influenced by a balance between matrix metalloproteinases (MMPs) and tissue 
inhibitors of metalloproteinases (TIMPs). A significant increase in the TIM P-1 level 
at 1 week after implantation compared with that in gingival crevicular fluid (GCF) 
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from healthy periodontium has been reported. Four weeks after implantation it had 
reached the same level as that in the GCF of healthy subjects. The data has also 
disclosed a higher post-implantation collagenase activity level at 1 week than those 
at weeks 2, 4, and 12. This may be due to the increase in MMP-1 and -8. The 
overproduction of TIMP-1 in the wound area after implantation could, to some 
extent, inhibit excessive tissue destruction and degradation of the neo-matrix in 
wound repair due to MMPs *. 

Four biomechanical parameters (peak force, vertical displacement, interface 
stiffness, and strain energy) were measured at implant placement and after 3 months 
of healing during push-in tests on implants placed in the mandibles of dogs. Results 
showed increasing values of peak force, interface stitfness, and strain energy after 3 
months of healing, but decreasing value of vertical displacement. These findings 
suggest that the interface stiffness, which is considered as a major factor for implant 
success, increases during 3 months of healing in dogs (correspondent to a 4 to 6 
month healing period in human mandibles), which indicates that a 4 to 6 month 
healing period for dental implants is justifiable.” 

Tissue damage due to madequate surgical techniques, resulting in excessive bone 
necrosis; bacterial contamination of the implant and/or the wound in course of 
surgery or insufficient maintenance care procedures in the healing phase of I-stage 
implants; and the instability of the implant and/or premature loading, are considered 
to be the key factors that can influence bone healing around titanium implants”. 
Studies on the effect of bacterial infection on the periimplant tissues revealed a 
decrease in woven bone formation due to infection”’. 


2.3 The Soft Tissue Interface 


Epithelial cells have the capacity to proliferate and to move on surfaces. Following 
implant/abutment placement, epithelial migration from adjacent soft tissues begins. 
The epithelium moves in apico-coronal direction as soon as it reaches the implant 
surface, giving rise to a junctional epithelium about 2mm long.” The formation of a 
junctional epithelium in the implant/mucosal interface can be considered the first 
barrier of defense against oral microflora. Research has shown that once epithelial 
cells have reached the implant surface, their attachment occurs directly via a basal 
lamina (< 200nm) and hemidesmosomes.”*’ Another form of attachment, where a 
glycoprotein layer of approximately 200nm forms exists between the cell wall 
without direct epithelium/implant contact was also proposed.*” Nonetheless, it is 
agreed that this seal is probably viable and adequate in function since there were 
minimal histologic inflammatory reactions in the underlying connective tissues. “”' 
The surface topography has also played a role in the soft tissue attachment to 
titanium surfaces, since polished surfaces have shown higher compatibility to 
fibroblasts when compared to rough surfaces.”*’ A prerequisite to a successful 
dental implant should obtain a peri-implant mucosal seal on the implant surface.*” 
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Failure to achieve or maintain this seal results in apical migration of epithelium into 
the bone/implant interface and possibly caused partially or completely encapsulation of 
endosseous implant.”°' In a natural dentition, the junctional epithelium provides a 
seal at the base of the sulcus against the penetration of chemical and bacterial 
substances. If the seal is disrupted and/or the fibers apical to the epithelium are 
destroyed, epithelia cells migrate apically, forming a periodontal pocket after 
cleavage of the soft tissue from the radicular surface. The importance of the biologic 
seal should be emphasized due to the absence of cementum and fiber insertion into 
the implant surface which may cause formation of a “ peri-implant pocket” 
extending into the osseous structures. 

For the connective tissue attachment to the implant, collagen fibers form a tight 
cuff around the implant abutment. The length of the connective tissue attachment 
ranges from 1.3 to 1.8mm””’, and is also dependent on the design of the implant 
(one- or two-part implant). There are reports in the literature that the use of a 
plasma sprayed surface promotes connective tissue adherence with fibers inserted 
functionally at 90° into the plasma sprayed surface of the implant.*” However, 
most reports have shown that connective tissue fibers are not inserted into the 
implant surface. Dense bundles of thick collagen fibers are oriented only 
longitudinally due to a lack of fiber insertion.*” The connective tissue of the peri- 
implant mucosa can be subdivided into two distinct zones. The outer zone located 
under the junctional epithelium, composed of collagen Types I and III, is 
responsible for the transformation of collagen. The inner supracrestal connective 
tissue zone, composed mainly of Type I collagen, is responsible for mechanical 
resistance and stability of the periimplant mucosa.” The connective tissue 
surrounding dental implants possess fewer capillaries as compared to that 
surrounding natural teeth.’ The blood supply of the peri-implant mucosa is 
composed of terminal branches of larger vessels originating from the periosteum of 
alveolar bone, while the vasculature of the gingiva and the supracrestal connective 
tissue at teeth is derived from the supraperiosteal vessels lateral of the alveolar 
process and the vessels of the periodontal ligament.*’ The peri-implant mucosa has 
also shown to exhibit an inflammatory response, similar to gingival tissues, 
characterized by increased proportions of T- and B-cells in the inner connective 
tissue. However, the host response is less pronounced than in gingival tissues. *® 

Animal and human studies have shown that an adequate biological width/seal can 
only be achieved if there is a supracrestal smooth titanium surface of at least 3 mm 
long in the apico-coronal direction.” The junctional epithelium(JE) occupies 
approximately 2 mm of this surface, while the rest of it is occupied by the 
connective tissue (CT). The values for two-part implants were 2.0 mm for JE, and 1. 
3 to 1.8 mm for CT.” The values for a biological width/seal around one-part 
implants at 3-month unloaded implants were reported to be 0.49 mm for sulcus 
depth (SD), 1.16 mm for JE attachment, and 1.36 mm for CT attachment. These 
dimensions slightly varied in the loaded implants, 0.50 mm for SD, 1.44 mm for JE, 
and 1.01 mm for CT. After a 12-month loading period, these values were 0.16 mm 
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for SD, 1.88mm for JE, and 1.05mm for CT.” It can be speculated that the biological 
width/seal increases approximately Imm after implant loading possibly due to 
crestal bone resorption. Table 2.2 summarizes the comparison between dental 
implants and natural teeth. 


Table 2, 2 Dental Implants vs. Natural Teeth 


Dental Implants Natural Teeth 
Hard Tissue Interface Only bone (functional ankylosis) Bone, PDL, cementum 
Hemidesmosomes and basal Hemidesmosomes, and 
Soft Tissue Interface ; 
lamina basal lamina 
CT Fiber Apparatus Parallel Perpendicular 
Collagen Fibers More Less 


2.4 Protocols for Implant Placement 


A dental implant is said to be biomechanically osseointegrated if no progressive 
relative motion of living bone and implant occur under functional levels and types of 
loading for the entire life of the patient. 

The surgical phase of most dental implant systems involves two distinctive 
stages, the placement of the implant and its subsequent uncovering. Research 
reported no significant differences in wound healing of implants exist in one-stage 
and two-stage procedures.” However, when two-stage implants were inserted 
using a one-stage technique, failure rate of 16. 7% was observed, suggesting that 
two-stage implants can be used in a one-stage approach but with increased risks of 
failure.” Primary stability of dental implant is considered to be a key factor to 
promote close bone-to-implant contact. Less bone to titanium contact and lower 
bone level in the threads have been determined in initially mobile implants compared 
to initially stable ones.” 


2.5 Types of Implant Systems 


Advances in oral implantology research have led to the development of several 
different types of implants, and it is anticipated that continued research will lead to 
even improved systems. The different implant systems are classified according to 
their morphology and method of osseous incorporation that are subperiosteal, 
transosteal, and endosseous implants. The most frequently used implants are 
endosseous implants while other forms have been slowly replaced by consumers. 
Endosseous implant systems include a range of sizes, shapes, coatings, and 
prosthetic components (Table 2.3). A variety of lengths and widths should be 
available to better incorporate the implant fixture within osseous structures. 
Prosthetic components can also be selected in a variety of sizes and angles to 
perfectly accommodate the final restoration. Implant surface morphology has been 
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shown to influence osseointegration. Porous coatings (i.e. acid-etched, sand-blasted) 
can achieve more bone-to-implant contact than smooth subcrestal surfaces.” Chemical 
coating of implant surfaces (i.e. hydroxyapatite) has also been shown to promote 
enhanced bone-to-implant contact and faster osseointegration. 


Table 2,3 Classification of Implant Fixtures According to Surgical Protocol, 
Interface Connection, and Body Geometry 


Surgical Protocol | Interface Connection Body Geometry 
~— One-stage — Internal connection — Threaded screws 
— External connection — Straight 

-— Tworstage — Octagonal — Tapered 
— Hexagonal — Stepped 
—— Cone Screw — Expansive 
=o 4 One T16X 
—~ Spline 


With an immense variety of available implant systems, the clinician should 
follow appropriate criteria for selection of an implant sy stem (Table 2.4). However, 
these systems appeared to be similar and very few have been subjected to a variety 
of clinical trials to determine their success. The majority of implant systems 
currently available were modifications of those few systems that have been 
properly evaluated. However, the data of one system should not be applied to 
another. It is important to know the longterm stability of the fixtures and the 
degree of tolerance of the peri-imp lant tissues for extended periods of time. An ideal 
implant system should also be easy-to-use and flexible to a variety of clinical 
scenarios. Only systems that provide numerous restorative options can predictably 
result in pleasing esthetic outcomes. 


Table 2. 4 Dental Implant Selection Criteria 
Selection Criteria 


— High success rate 

— Predictable Osseointegration demonstrated in scientific data 
— Prosthetic flexibility in a variety of clinical applications 

— Tissue friendly 

— Interface stability 

— User friendly 

-— Optimal esthetics 

— Reasonable cost 


2.6 Prosthetic Rehabilitation 


The exact position and angulation of the implants in jaw bone is planned before 
surgery. Utilization of surgical guides are often needed to precisely locate ideal 
implant position during placement. Varying degrees of edentulism can be treated by 
prosthetic rehabilitation with implants. These range from a fully edentulous jaw to 


50 Application of the Finite Element Method in Implant Dentistry 


partially edentulous segments to single tooth replacements.(Fig.2.2a-2.2)). 
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Fig. 2.2b Radiograph Demonstrates Adequate Bone Height for Implant Placement 
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Fig. 2.2¢ Implant Osteotomy Is Prepared under Guidance of a Surgical Stent 
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Fig. 2. 2d Intra-surgical Radiograph Shows Ideal Angulation for Implantation 


Fig. 2. 2e Implant Is Placed in a One-stage Approach with Satisfactory Position 


Fig. 2. 2f Interrupted Sutures Reposition the Flaps 
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Fig. 2. 2g Post-surgical Radiograph Confirms Adequate Relation 
between Implant and Adjacent Teeth 


Fig, 2.2h Adequate Healing Is Observed 2 Weeks Later 
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Fig. 2. 2i Four Months Later, the Final Crown Is Cemented 
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Fig. 2.2j Radiograph Demonstrates the Final Restoration 
Complete dentures are always confronted with the problems of retention and 
stability, especially in patients with advanced alveolar ridge resorption. Implant 
supported dentures allow for much improved masticatory function and comfort. 
Commonly, the lower complete denture is the more problematic and thus it 
frequently indicates for implant supported rehabilitation. (Fig, 2.3a — 2.3e) 


Fig. 2. 3a Radiograph Demonstrates Edentulous Maxilla with Insufficient 
Alveolar Bone Height for Implant Placement 


Fig. 2. 3b Following Augmentation of the Maxillary Sinuses, Four Implants Are Placed 
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Fig. 2.3e Satisfactory Esthetics and Occlusal Function Are Achieved 
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Even when the patient retains some teeth, some edentulous spans are too long 
for predictable treatment with conventional bridges. When the edentulous span has 
no terminal abutment tooth at the end of the arch to support a bridge, a so-called 
free-end saddle situation is said to exist. Such patients had no choice but to wear 
removable partial dentures in the past. Implant supported bridges are now possible 
in these situations. Increasing numbers of patients are opting for the restoration of a 
single or multiple adjacent missing teeth with implant-supported crowns. In 
conventional dental treatment, the teeth adjacent to the edentulous space have to be 
altered to serve as abutments for fixed partial prostheses. If these abutment teeth 
are untouched by previous episodes of caries then such a treatment modality would 
be destructive to healthy tooth structure. Before the advent of implants, this was 
the best restorative option for partial edentulism. Replacement of these teeth with 
imp lant-supported crown has become the standard of treatment in the last decade 
(Table 2.5). 


Table 2.5 Dental Implants vs. Natural Teeth 


Dental Implants Natural Teeth 
Hard Tissue Interface Only bone (functional ankylosis) Bone, PDL, cementum 
Sch Tiseus Tntahce Hemidesmosomes, and Basal Hemidesmosomes, and 
Lamina Basal Lamina 
CT Fiber Apparatus Parallel Perpendicular 
Collagen Fibers More Less 
Soft Tissue Vascularity Less More 
Fibroblasts Less More 
Probing Depth in Health Deeper (3-4 mm) Shallower (< 3mm) 
Inflammation Progression Faster Slower 
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5.1 Introduction 


As the Finite Element (FE) method requires a huge amount of computation, its 
application should be supported by advanced of computer technology. Many 
computer softwares have been developed and commercialised for workstation and 
Personal Computers (PC) in recent years. Most of the FE softwares are developed 
for general applications in solids, fluids, and electrical and magnetic field. There are 
also some FE software designed for a specific type of analysis. The FE software 
used for dental implant simulation is the one for general applications. 

The selection of FE software is based on requirement of applications. Although 
the fundamental principle of FE analysis is the same, computer softwares can be 
developed with different features. Some commercial FE softwares have well 
developed CAD function, while others provide more problem solving options. A 
simple dental implant FE simulation is typical solid analysis, which can be 
performed using most of currently available commercial FE softwares. However, 
when the osseointegrated interface takes into account, contact analysis is required. 

In this part, two well-known FE softwares, ANSYS and ABAQUS, are used as 
examples to illustrate how the dental implant simulation is performed using a 
commercial FE software. Users who are new to the FE analysis could become 
familiar with the software through the step by step exercise described in the 
following chapters. 
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5.2 ANSYS 


5.2.1 Introduction 


ANSYS is one of the most popular commercial FE softwares. It is widely used to 
solve problems ranging from typical mechanical engineering to medical engineering in 
engineering and academic fields . It can be adapted to simulate simple linear analyses 
as well as the most challenging nonlinear dynamic analysis. ANSYS contains an 
extensive class of element types, which could virtually model any geometry and an 
equally extensive list of material properties. In terms of mechanical analysis, it not 
only can simulate the behaviour of most typical engineering materials including 
metals, rubber, polymers, composites, crushable/resilient foams, etc, but also 
biomaterials such as soft tissue and bones. 

ANSYS is a user-friendly software and very simple to use. Compared to other 
commercial FE softwares such as ABAQUS, it not only offers users a wide range of 
capabilities but also powerful 3D modelling functionality. For example, users can 
easily create an engineering geometry within ANSYS. This is particularly useful for 
users without engineering background. In addition to its command line interface, 
ANSYS provides a Graphic User Interface (GUI). The GUI is featured in a style of 
step-by-step menu. Following this well structured menu, a FE model can be created 
by providing the geometry data, material properties, boundary conditions, and the 
loads applied to the model. Once the model is created, to solve the simulation 1s 
straightforward in ANSYS. Most importantly, ANSYS has very effective post 
process capabilities from displaying colour strain contour to plotting stress 
distribution on a user-defined path. Users rarely have to define parameters for 
controlling the numerical solution of a problem. 

ANSYS is modularised a “standalone” package with three fundamental modules. 
They are Preprocessor, Solution and General Postprocessor modules. These modules 
are under the ANSYS Main Menu of the ANSYS GUI. In addition, there are some 
specialised modules such as Time History General Postprocessor, Topological 
Option and Radiation Option. FE implant simulations are covered by the first three 
fundamental modules. Therefore, In this chapter, we will concentrate on these 
fundamental modules. 


Preprocessor 


ANSYS Preprocessor module is the one used for creating a FE model. It consists of 
a number of sub-menus, which can be used for geometrical model creation, assigning 
element type and material properties, mesh generation, and applying constrains and 
loads. In the case of a complicated model geometry, ANSYS can import the model 
created by other softwares, including IGES, STEP, Pro Engineering, and Parasolid. 
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Solution 


ANSYS Solution module is the one used to solve the FE model. There is some 
overlap between this module and Preprocessor. Either of them can apply loads and 
constrains to the model. The main function of the Solution module allows users to 
configure how they want the model to be solved, like the Analysis Type users used, 
how the load steps are applied to the model. 


General Post processor 


ANSYS General Postprocessor module is the one used to review the solution result 
of the FE model. It provides some very useful tools to view results of the solution. 
It includes plotting and listing stress, strain and displacement contours of the model. 
With plot on graphic function of the module, users can define a path on the model 
and plot a result (stress, strain or displacement values) on the path. Both graphic 
and numerical data can be easily exported and saved into a format that can be read 
by other softwares such as MS Excel for further process. 

In order to provide a guideline for users who have no engineering background to 
learn how to use ANSYS for dental implant simulation, a step-by-step FE modelling 
example of a dental implant is described in the following sections. 


5.2.2 Preprocess 


The creation of a FE model is done by Preprocessor module. Once the geometry of 
the model is defined (in most cases, the geometry comes from a real application), the 
first thing is to convert it into a FE model. Assumptions need to be made, as the FE 
model is never to be the exact same as the real case. This may require some 
experiences and practises. The assumptions given by this step-by-step case study 
can be taken as a starting point for creating a FE model of dental implant. 

The case used in this chapter is to simulate a biting force applied on a 
Branemark “same day” dental implant. The aim of this study is to investigate the 
stress distribution in the bone and implant. To simplify the study, simple bone and 
implant structures are used, and following assumptions are made: 

(1) good osseointegration exists between bone and implant, i.e. no relative movement 
between bone and implant occurs under loading condition; 

(2) the implant components are treated as a continued structure; 

(3) bone is treated as homogeneous material in terms of elastic module; 

(4) the mandible is fully constrainted on its ends. 


5.2.2.1 Selection of Element Type 


The first step of FE simulations is to select element type. The selection depends on 
the case of study. To add element type to the model, click Add/Edit/Delete under 
the Element Type menu, a Element Type Window will pop up as shown in Fig5.1. 

The Add button in the window shown in the Fig.5.1 is for creating new element 
type to the model. Click Add button, a Library of Element Types window pops up 
as shown in the Fig.S.2. 
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Element Types 


Defined Element Types: 


NONE DEFINED 


MB Library of Element Types 


Library of Element Types Structural Mass A |Brick Snode 45 

8node 185 
20node 186 
20node 95 
layered 46 
layered 191 

Hyperelastic 

Mooney-Riviin 


Element type reference number 


Fig. 5.2 Library of Element Types Window 


For a stress analysis model, a solid structure element type should be used. As 
shown in the Fig.5.2, there are more than one options of the solid structure element 
type such as Brick 8node 45, 8node 185. Tetro Hedro element type can be used for 
structure analysis as well. The general guideline for selecting element type is: for 
regular shape such as a block, a cylinder, or a uniform section beam, Brick type 
element is normally used; for irregular geometry, such as 3D model of mandible, 
Tetro Hedro element type is used, as it can fit to the geometry well. In this 
example, a Tetro Hedro element type-10node Element 92 is used. To select this 
element type, it needs to drag the scroll bar down to the bottom of the window to 
see the element type. Then use the mouse to highlight the element before clicking 
the OK button. After the OK button is clicked, a Element Types window pops 
back and SOLID 92 element appears in the window as shown in the Fig.5.3. 

To add more element types, clicking add button again. To quit the Element Type 
window, clicking the Close button. 
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Element Types 


Defined Element Types: 
Type 1 SOLID92 


Fig. 5.3 Element Type Window with Inserted Element Type SOLID92 


5.2.2.2 Assign the Material Properties 


The second step of ANSYS FE simulation is to assign material properties to the 
model. It can be done before or after the creation of model geometry. For stress/ 
strain analysis, there are two essential parameters need to be defined, Elastic 
modulus and Poisson’s Ratio. To define a material property, click the Material 
Models under Material Properties menu, then a Define Material Model Behavior 
Window pops up as shown in Fig5.4. 


A Define Material Model Behavior 
Material Edit Help 


Material Models Defined Material Models Available 


© Material Model Number 1 o Structural 


$y Thermal Expansion Coef 
@ Damping 
© Friction Coefficient 


Fig. 5.4 A Define Material Model Behavior Window 


In the Define Material Model Behavior Window, users can define linear, 
nonlinear, isotropic, orthotropic, and anisotropic properties to the model according 
the nature of the simulations. In this example, the bone is assumed to be isotropic 
and loading is in a linear range, i.e. to define linear isotropic properties to the model. 
Click the Isotropic, a new Linear Isotropic Properties for Materials Number 
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window pops up as shown in Fig.5.5. In this window, users can fill Elastic modulus 
and Poisson's Ratio. In this example, it is assumed that Elastic modulus and 
Poisson's Ratio of cortical bone are 18GPa and 0.3, respectively. 


a ee 


Linear Isotropic Properties for Material Number 1. 


Fig. 5.5 A Linear Isotropic Properties for Materials Number Window 


To add another material property, click an appropriate type of material 
behaviours in the Define Material Model Behavior window again. After all material 
properties are defined, users can close the Linear Isotropic Properties for Materials 
Number window by click OK and Define Material Model Behavior window by click 
the small red cross on right top corner. 


5.2.2.3 Model Creation 


The third step of the FE simulation is to create model geometry. In ANSYS, a 
simulation can be carried out on 2D or 3D geometry. A component with simple 
geometry can be simplified as a 2D case, such as a uniform cross section beam. 
Single dental implant simulation can also be carried out on a 2D geometry and 
achieve good approximation. The guideline of selection the model geometry is unless 
it is necessary to use 2D simulation. It will save both modelling and computation 
time of the simulation. 

To study the same-day-teeth implant system, it deals with bridge structure of 
dental fixture and can not be modelled using 2D geometry. Therefore, a 3D model 
has to be used. As an example, an approximate 3D model is created in this section to 
show users how to create a 3D model in ANSYS. This model is far more close to a 
real same-day-teeth implant system. 

To create a FE model, two critical questions need to answer, how accurate does 
a model need to be and where does the geometry information come from? These are 
not straightforward questions. To answer them, users should firstly identify what 
they want to get from the simulation. Is the local stress around the dental implant, 
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deformation of the bridge fixture, overall stress in the bone or the stress in the 
implant? In this example, the deformation of the bridge fixture and approximate local 
stress around the implant are what we are interested in. Therefore, a simplified bone 
structure together with an approximate geometry of the dental implant system is 
used. 

It is very obvious that the more accurate the model is, the more difficult it will 
be in creating a model. For approximation and simplification, in this example, the 
dental implant is modelled without thread. In a more accurate model, the implant 
needs thread and the bone geometry information should come from human anatomy, 
i.e. CT scan images, which is shown in the next example of this chapter. 

For model generation using human anatomy data, ANSYS can import solid model 
in different formats, such as IGES, STEP, Pro Engineering, Parasolid and so on. 
ANSYS also can read text file written in its own language, by doing this it can 
import coordinates from CT images, from which a more accurate model is created. 

This example concentrates on how to create a simple dental implant-bone model 
as shown in Fig.5.6. The model shows a bridge fixture supported by three 3mm 
implants. The length of the implant is 12mm, the thickness of bridge plate is 3.5mm, 
and the depth of the bridge plate is 6mm. In the model, the mandible is modelled as 
a single piece with thickness of 6.5mm. 


Fig. 5.6 A Simplified Dental Implant Bone Model 


Creation of 3D model in ANSYS 

Although there are special commands can be used to create simple geometry such as 
cylinder and brick, to create a 3D model in ANSYS, such as bridge plate, following 
steps are still recommended: create Key points from the coordinates of the plate 
geometries, create Line from the Key points, create Area from the Lines, and finally 
create volume from the Areas. 
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Key point is defined by coordinates in 
ANSYS, which is the lowest level of 
geometry. Unless ANSYS provides a 
special command to directly create a 2D or 
3D geometry, model generation starts from 
creating Key points and follows the order 
described above. To create a Key point, 
click the Create command under the 
Modeling menu. A list of geometry 
creation command will appear including, 
Lines, Areas, Volumes, Nodes, and 
Elements. The Keypoints command is the 
first one as shown in Fig.5.7. 

Under the Keypoints command there 
are 7 options to create Key point in 
ANSYS. Among these options, the In 
Active CS (Coordinate System) is 
commonly used when the coordinates of 
the model are provided. 

Click the In Active CS in the menu, a 
Create Keypoints in Active Coordinate 
System window appears as shown in 
Fig, 5. 8. This window provides space to 
enter the Key point number and 
coordinate values. The Key point number 
is the identification of each Key point. 


AAA ANSYS MainMenu AAA 


sferences + 
eprocessor 
Element Type 
Real Constants 
’ Material Props 
Sections 
Modeling 
E)] Create 
= Keypoints 
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Fig. 3.7 Keypoints Creation Menu 


Each Key point has an unique Key point number. 


@® Create Keypoints in Active Coordinate System 
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Fig. 5.8 A Create Keypoints in Active Coordinate System Window 
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Please note that the ANSYS uses ISO unit. The value appears in the coordinate 
textboxes are in meter. ANSYS has a command Scale can be used to re-scale the 
model. This is very useful especially when a model is imported from other CAD 
software or format of 3D solid model, where meter is often used as the unit. The 
Scale command is used to reduce the model by 1/1000. 
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Once necessary numbers of Key point are 
created, the next step is to join the Key 
points by Lines in the way of model goes. 
There are three types of lines: Straight line, 
Arc and Spline. A straight line is defined by 
two Key points, a Arc by three Key points, 
and a Spline by more than three Key points. 
Users can select them according to the 
requirement of the model. 

To create a straight line, arc or spline 
through Key points, click the Lines command 
under the Create menu. A list of Line 
command will appear in the window as 
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shown in Fig.5.9. 


There are other commands in this menu 
which allows users to create lines. 


Splines 


Fig. 5.9 Line Command Menu 
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Fig. 5.10 Create Straight Line Window 


Straight Line command, a Create 
Straight Line window pops up as 
shown in Fig5.10. Pick two Key 
points and join them using a straight 
line as shown in Fig.5.11. 

In a simple model, users can create 
all Key points before join them by 
lines. But in a complicated model, it is 
recommended not to create all Key 
points at a time, otherwise it will be 
too confused to pick Key points to 
form lines. Users need to have a good 
plan in creating a complicated model. 
Experience is very important to master 
the skill of model creation in ANSYS. 

After lines are created, the next 
step is to create areas from lines. To 
create an area, click the Areas 
command in ANSYS Main Menu 
shown in Fig.5.7. Under the Areas, 
there is an Arbitrary command. In 
ANSYS, area can be created by both 
Key points and lines. It is 
recommended to use lines to create 
areas because this would generate less 
confusion. 
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Fig. 5.11 A Straight Line Join Two Key Points 


Click the By Lines command under the Arbitrary command, a Create Areas by 
Lines window pops up.Pick option is selected by default. Users can pick the lines 
in the model window to form an area as shown in Fig.5.12. 


Fig. 5.12 Creating Areas by Lines 


An area geometry can be used for 2D FE meshing and modelling. For 3D FE 
simulations, although a model can be meshed through Key points of which the model 
consists, meshing through volumes are much easier than using Key points of which 
the volume consists. And some Boolean operations such as volumetric Intersect, 
Add, Subtract, and Glue algorithm can only be carried out between volumes. 
Therefore, it is strongly recommended for beginners to follow the steps given here 
to create FE models. 

To create a volume by areas, click the Volume command in ANSYS Main Menu 
shown in Fig.5.7. Under the Volume, there are six volumetric creation commands. 
Only the Arbitrary command allows users to create volumes from areas. A volume 
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can also be created by Key points. It is obvious, however, to create it by areas is 
much easier than directly by Key points. Fig.5.13 shows a simplified volumetric 


model of the same-day-teeth implant system. 


Fig. 5.13 A Simplified Volumetric Model of the Same-day-teeth Implant System 


5.2.2.4 Mesh Generation 


Mesh generation is the forth step of FE simulation. In this step, a 2D or 3D model 
is meshed with elements defined in the first step and material properties defined in 
the second step. The mesh process is to divide the geometrical model created in the 
third step into small and finite divisions. An element consists of a number of edges 
and nodes, which depends on the element type selected. During the mesh generation, 
material properties of the model are assigned. It 1s important to make sure each part 
of the model has been assigned right element type and material property. 


To assign element type and material 
property to a model can use Default Attribs 
command as shown in Fig.5.14 before meshing 
the model. 

The element type and material property 
can also be assigned to volume using Picked 
Volumes command. Click the command, a 
Volume Attnbutes window pops up. Click the 
volumes in the model window. After click the 
OK button a new one as the same pops up as 
shown in Fig.5.15. 

In this window, users can select material 
number and element type number to be signed 
to the volume picked up. 

There are two ways to generate mesh in 
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Fig. 5.14 Meshing Menu 
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Fig. 5.15 Volume Attributes Window 
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adapted with fine element at these places. But in 

Other situations that the stress may not vary much, a coarse mesh would be 
sufficient enough. In ANSYS, these can be achieved using Smart Size, as it takes 
model geometry into account while the mesh is generating. 

For manual mesh, users need to define element length on the geometry. This can 
be done by defining the number of elements or actual size on lines. To achieve an 
adequate mesh generation, experience is essential in this situation. 

The beginner of ANSYS is recommended to use Mesh Tool and select Smart 
Size. To use Mesh Tool, click the MeshToo! command in the menu as shown in Fig. 
5.16, a MeshTool window pops up as shown in Fig,5.17. 

In the MeshTool window, users can select mesh shape, free or mapped mesh, 
and option of Smart Size. As it shown in Fig.5.17, the size of element can be 
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© Fea 


Fig. 5.17 A MeshTool Window 


Mesh | 


adjusted between 1 and 10 when the check box of the Smart Size 1s ticked. After all 
the settings are selected, click the Mesh button to start meshing, and then a Mesh 
Volumes window pops up. Pick the volumes to be meshed in the model window. A 


meshed model is shown in Fig.5.18. 
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Fig. 5.18 A Meshed Implant System 


5.2.2.5 Loads and Constraints 


Once the model is meshed, the next step is to apply structural loads and constraints 
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to the model. In dental implant FEA, the structural loads are forces or pressures 
being put on the implant or bridge fixture and the constraints are boundary 
conditions of how, for example, the mandible is restrainecl under loading. 

The structural loads can be defined in ANSYS using point force, pressure on one 
area, and pure moment. However, not all element types allow users to apply a pure 
moment to the model. The element selected for this dental implant analysis is 
SOLID92. Only point force and pressure can be applied to the FE model. 

To define the boundary conditions of the model, we assumes that the movement 
of the mandible is fully restrained in all directions under loading and the constraint 
is located at two ends of the mandible as shown in blue colour in Fig.5.19. 


Fig. 5.19 The Restrain Locations on the Mandible 


To apply loads on the model, users select force or pressure according to the 
requirement of simulation. The force is applied to Nodes, Keypoints or Areas. The 
pressure 1s only applied to an area. In ANSYS, loading and constraints can be done 
by either Preprocessor or Solution module. | 

To apply a force load, click Define Loads under the Load of the Preprocessor 
module as shown in Fig. 5.20 (a). In the ANSYS Main Menu, click Apply command 
to apply structural loads and constrains to the model. Once the Apply is clicked, 
more options appear in the ANSYS Main Menu as shown in Fig. 5.20(b). The dental 
implant FE simulation is a structural analysis, click the Structural command to bring 
more options as shown in Fig. 5.20 (c). 

To apply constraints to the mandible, click the Displacement command, more 
options appear. These options allow users to apply displacement constraints on 
Nodes, Keypoints, Lines and Areas. Click the Areas command in the main menu, an 
Apply U,ROT on Areas window pops up as shown in Fig,5.21(a). 

Click the areas on the model to highlight the areas to Apply U,ROT on Areas 
window be constrainted. Once all areas are selected, click the OK button in the 
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Fig. 5.20 ANSYS Main Menu for Applying Loading 
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Fig. 5.21 Apply U,ROT on Areas Windows 


Apply U,ROT on Areas window, a new Apply U,ROT on Areas window as shown 
in Fig.5.21(b). In this window, users can define the displacement in three directions, 
individually or together. In the dental implant model, no movement has been 
assumed in all three directions. To define this, select All DOF (Degree of Freedom) 
in the window and enter 0 in Displacement value box, and then click OK button to 
complete applying constraints to the model. 

To apply a load (force or moment) to the model, click the Force/Moment or 
Pressure command. As shown in Fig. 5.22 (a), there are more options for applying 
load. The force/‘moment can be applied to Keypoints or Nodes. In ANSYS, the 
moment can only be applied to certain elements with rotational DOF. The pressure 
can be applied to Lines, Areas, Nodes or Elements. 

To apply force on the dental implant, click On Nodes command in the ANSYS 
Main Menu, an Apply Force/Moment on Nodes window pops up as shown in Fig. 
5.22 (b). Use the mouse to select one or more nodes on the implant model. Click the 
OK button in the Apply Force/Moment on Nodes window, a new Apply Force/ 
Moment on Nodes window pops up. In this window, users can select which 
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Fig. 5.22 Menu for Apply Load to ANSYS Model 


direction the force will be applied. Click the arrow button, a dropdown list appears 
as shown in Fig.5.23. In this list, users can select one direction from x, y and z. Once 
the direction is selected, enter a value in the Force/moment value box. At this point, 
if users click OK button, the load application completed. If the Apply button is 
clicked, users can continue to apply force or moment to the model. To delete the 
applied load, click the Delete / Structure / Force/moment command. 


Fig. 5.23 Dropdown Window for Apply Load in Different Directions 


As shown in Fig.5.24, Once the load is applied, a red arrow indicator will appear 
at the point where the force is applied. 


Fig. 5,24 Load Indicator on Implant Model 


5.2.3 Solution 


Once the mesh is completed and boundary conditions are applied, the next step is to 
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solve the model using solution module. This module starts Analysis type. There are 
three options under Analysis type, New Analysis, Restart, and Solution control. 
For the first time analysis, use the New analysis. After click the New analysis 
command in the ANSYS Main menu, a New Analysis window pops up, as shown in 
Fig.5.25. In this window, users can select the Type of analysis. For the dental 
implant simulation, the Static analysis type is selected. 


MB New Analysis 


[ANTYPE] Type of analysis 
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(~ Spectrum 
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€~ Substructuring 


OK | Cancel | Help | 


Fig. 5.25 New Analysis Window 


The Solution Control option allows users to assign different parameters to the 
model. For static linear analysis, the default settings can be used. For other analysis 
such as nor linear and transient, users need to properly define these parameters. For 
further information about how to define these settings refers to ANSYS manual. 

Before the model is solved, users also need to specify the loading steps and 
output format of the solution. This is done using Load Step Option command. For 
multi-step analysis, users need to define Printout format of the result. Click the 
Solution Control command in the ANSYS Main Menu, a Solution Printout Controls 
window appears as shown in Fig.5.26. Users can define the solution printout using 
the options provided in this window. 

The final step of the model solution is to execute a sub-command under the 
Solve command. The Current LS command is one used for the first time analysis. It 
assigns the current load setting to the model, and then carry out the computations to 
solve the model. 


5.2.4 Postprocess 


Once the model is solved, results of the analysis can accessed using General 
Postprocess module in ANSYS. For structure analysis such as a dental implant 
system, this module provides users with three fundamental functions to review the 
results. They are Plot Result, List and Export Result, and Plot graphs. 


5.2.4, 1 Plot Result 


The Pilot result function allows users to generally review the results of analysis in a 
format of contour or vector graph. For structural analysis such as a dental implant, 
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Fig. 5.27 Contour Nodal Solution Data 


the stress distribution, strain level, and deformation of the structure can be shown 
using contour or vector plot. 

To plot stress distribution contour, click the Node Solution under the Contour 
Plot command in the ANSYS Main Menu. A Contour Nodal Solution Data window 
pops up, shown in Fig. 5.27. 

In the Contour Nodal Solution Data window, users can select parameters to be 
plotted such as stress, strain or displacement. There are many options under each of 
the above parameters. Fig.5.28 (a) and (b) shows that under the stress parameter, 
there are options of stress in x , y and z directions, shear stress in xy, yz, anc xz 
directions, principal stresses, von Mises stress, and so on. 
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Fig. 5.28 Sub-option of Stress Contour Plot 


To plot von Mises stress contour, use the mouse to highlight the von Mises 
ESQV option shown in Fig. 5.28 (b). Once all option m the Contour Nodal Solution 
Data window has been defined, click the OK button in the window, a von Mises 
stress contour will appear in the window as shown in Fig, 5.29. To plot strain and 
displacement contours is similar as to plot stress contour. The procedure of vector 
plot is similar as the contour plot. The only difference between the vector and 
contour plots is that the former uses vector indicators (small arrows) to represent 
the results, while the later uses colours to show the result. 

The advantage of contour format is to give quick and direct graphic view of the 
result. Its limitation is that it only shows the result on the surface of the model and 
does not provide any accurate value and location. 


5.2.4.2 Plot Result on Graph 


Plot result on graph is a useful tool in ANSYS. It allows users to define a path 
across the model and plot the analysis result against the defined path to a graph. By 
using it, users can specify the location on the model and review the result accurately 
at this specific location. To plot analysis result against a path across the model, 
users need to define a path first. Secondly, map a specific result to the path before 
finally plot the result on the path. In ANSYS, the Define Path command is under 
Path Operations command as shown in Fig. 5.30 (a). Path can be defined by nodes, 
on working plane, or by location. To define path by nodes, click the By Nodes 
command, and then a By Nodes window pops up (Fig. 5.30 (b)). Use mouse to 
select a node on the model as the start of the path. Users can select two or more 
nodes to define a path. Once all nodes are selected, click OK in the By Nodes 
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Fig. 5.29 von Mises Stress Distribution in the Dental Implant System 


window, the path define 1s completed and a new By Nodes window pops up as 
shown in Fig. 5,30 (c). 
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Fig. 5.30 Define Path by Nodes 


Enter a path name and set number of data or number of divisions on the path. 
The default values are 30 and 20, respectively. Click OK to confirm the setting in 
this window. The next step is to map the analysis result to the path. Click the Map 
onto Path command, a Map Result Items onto Path window pops up as shown in 
Fig. 5.31. 

In this window, users may enter User label for item and select the Item to be 
mapped as shown in Fig. 5.32. Once the settings in this window are defined, click 
OK to confirm the settings. 

To plot the mapped result on the path, click the On Graph command under the 
Plot Path Item, a Plot of Path Items on Graph window pops up as shown in Fig5.33. 
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Fig. 5.32 Define Settings in the Map Result Items onto Path Window 
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Fig. 5.33 A Plot of Path Items on Graph 


Select the entered User label for item (it is VMS in this example) and click OK, 
a graph will appears on the window as shown in Fig. 5.34. 


5.2.4.3 List and Export Result 


Although ANSYS provides users with powerful graphic plotting function, its list 


and export result function allows users to carry out process using spreadsheet 
software such as Excel. 


To list the result, click the Nodal Solution under the List Results command, a 
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Fig. 5.34 Plot of the Result Mapped to a Path 


List Nodal Solution window pops up as shown in Fig. 5.35. In this window, users 
selects the result to be listed. It is similar to selecting parameters to plot analysis 
result. Click OK after the parameter is selected; the result will appear in a popped 
up window (shown in Fig.5.36). 


MM List Nodal Solution 
[PRNSOL] List Nodal Solution 


Item,Comp Item to be listed i @ | |Components SCOMP 
Principds SPRIN 
Strain-total 
Nenlinear items 


Strain-plastic 
Strain-creep be PPrincipas SPRIN 


[AVPRIN] Eff NU for EQV strain 


Help | 


Fig. 5.35 List Nodal Solution Window 


To export the listed result to a file, click File, then Save as on the PRNSOL 
command window. This allows users to save a text file for later use. 


5.2.5 Summary 


In this chapter, a step by step description of how to use commercial finite element 
software ANSYS to carry out FE analysis for a Branemark same-day dental implant 
system is presented. It covers all the fundamental skills. In reality, the implant 
model may be more complicated than the one modelled in this chapter. However, the 
case shown in this chapter is a good starting point for those who want to use finite 
element analysis in their study. 


114 


6 
8 
8 
8 
6 
3) 
8 
8 


1 SUBSTEP= 1 


j088 


LOAD CASE= 8 


+ ARE FOR MATERIAL 1 


Si 


-14968E+87-6. 


-48716E +86 
- 93842E +86 
-19230E+87 
-21253E+8? 
-33392E+8? 


-51284E+86 


-8 


$2 $3 
66373E+86~-6.88193E+87 


-12048E+66-6.18176E+07 
6. 


16818E+86-6.49458E+67 


-27758E+86—-8.18192E+06 
-77226E+66 86261. 
-13523E+*8? S@715. 
-1442GE+6? 8.33819E+86 
-17886E+67-6.85674E+86? 


SINT 


6.180316E+68 
@.23847E+87 
6.58762E+8? 
8. 
4) 
4) 
4) 


26258E+07 


-20391E+67 
-32885E+8? 
-37162E+86?7 
8. 


96862E+8? 


Application of the Finite Element Methed in Implant Dentistry 


SEQU 
-94337E+8? 
-20748E+87 
-55113E+8? 
-19875E+87 
-17971E+87 
-28685E+8?7 
-3385SE+87 
-81763E+8? 


3138.9 
A 4AR ARE +A6 


-39877E +86 -8 .38326E+867 
1RNRIPIF+AP—-A RKRI11 F+AP 


8.38357E +8? 
A_99991 F +A? 


- 36547E+87 
RISVSE+AaAP 


4) 
8 
8 
-40543E+87 @ 
a) 
8 
—Aa 


Fig. 5.36 A List Results Window 


Comparing to other commercial software, ANSYS has demonstrated excellent 
suitability for various applications and easy to learn and use. Furthermore, the 
build-in preprocess and postprocess modules are very helpful for new users to learn 
how to use finite element software. 


5.3 ABAQUS 


5.3.1 Intreduction 


ABAQUS 1s a suite of powerful engineering simulation programs, based on the finite 
element method, that can solve problems ranging from relatively simple linear 
analyses to the most challenging nonlinear simulations. ABAQUS contains an 
extensive library of elements that can model virtually any geometry. It has an 
equally extensive list of material models that can simulate the behavior of most 
typical engineering materials including metals, rubber, polymers, composites, 
crushable/resilient foams, and so on. 

ABAQUS is a simple-to-use software with a wide range of capabilities, even the 
most complicated problems can be modeled easily. For example, problems with 
multiple components are modeled by associating the geometry to define each 
component with the appropriate material models. In most simulations, even highly 
nonlinear ones, users only need provide the engineering data such as the geometry of 
the structure, its material behavior, its boundary conditions, and the loads applied to 
it. In a nonlinear analysis, ABAQUS automatically chooses appropriate load 
increments and convergence tolerances. Not only does it choose the values for these 
parameters, it also continually adjusts them during the analysis to ensure that an 
accurate solution is obtained efficiently. Users rarely has to define parameters for 
controlling the numerical solution of the problem. 

ABAQUS consists of two main analysis modules—ABAQUS/Standard and 
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ABAQUS/Explicit. There are also two special-purpose add-on analysis products for 
ABAQUS/‘Standard— ABAQUS/Aqua and ABAQUS/Design. In addition, ABAQUS/ 
Safe provides fatigue postprocessing, while ABAQUS/ADAMS, ABAQUS/CAT, 
ABAQUS/C-MOLD, and ABAQUS/MOLDFLOW are interfaces to ADAMS/Flex, 
CATIA, C-MOLD, and MOLDFLOW, respectively. ABAQUS/CAE is the 
complete ABAQUS environment that includes capabilities in creating ABAQUS 
models, interactively submitting and monitoring ABAQUS jobs, and evaluating 
results. ABAQUS/Viewer is a subset of ABAQUS/CAE that includes just the 
postprocessing function. The relationship between these modules is shown in Fig. 
5.37. 
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ABAQUS/Standard 


ABAQUSStandard is a general-purpose analysis module that can solve a wide range 
of linear and nonlinear problems involving static, dynamic, thermal, and electrical 
response of components. This module is the main subject of this guide. 


ABAQUS/ Ex plicit 


ABAQUSExplicit is a special-purpose analysis module that uses an explicit 
dynamic finite element formulation. It is suitable for short-lived, transient dynamic 
events, such as impact and blast problems, and is also very efficient for highly 
nonlinear problems involving changing contact conditions such as forming 
simulations. ABAQUS/Explicit is not specifically covered in this guide but is 
described in the companion publication Getting Started with ABAQUS/Explicit. 
ABAQUS/CAE interface is the same for both analysis modules. In addition, the 
output is similar for the two modules, and ABAQUS/Viewer can be used to 
postprocess the results from either module. 


ABAQUS/CAE 
ABAQUS/CAE (Complete ABAQUS Environment) is an interactive, graphical 
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environment for ABAQUS. It allows models to be created quickly and easily by 
producing or importing the geometry of the structure to be analyzed and 
decomposing the geometry into meshable regions. Physical and material properties 
can be assigned to the geometry, together with loads and boundary conditions. 
ABAQUS/CAE contains very powerful options to mesh the geometry and to verify 
the resulting analysis model. Once the model is completed, ABAQUS/CAE can 
submit, monitor, and control the analysis jobs. The Visualization module can then be 
used to interpret the results. ABAQUS/CAE is also discussed in this guide. 


5.3, 2 Model an Implant in ABAQUS/CAE 


This example of a GJP biomechanical optimum implant, shown in Fig. 5.38, leads 
readers through the ABAQUS/CAE modelling process by visiting each of the 
modules and showing readers the basic steps used to create and analyze a simple 
model. 
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Fig. 5.38 ABAQUS/CAE Interface 


Part 

Sketch the two-dimensional geometry and create a part representing the implant. 
Property 

Define the material properties and section properties of the implant. 

Assembly 

Assemble the model. 

Step 

Configure the analysis procedure and output requests. 
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Load 

Apply loads and boundary conditions to the implant. 
Mesh 

Mesh the implant. 

Job 

Create a job and submit it for analysis. 
Visualization 

View the results of the analysis. 


5.3.2.1 Create a Part 


We use the Part module to create each of the parts to analyze. Parts that define the 
geometry of individual components of the model and, therefore, are the building 
blocks of an ABAQUS/CAE model. We can create parts that are native to 
ABAQUS/CAE, or import parts created by other applications either as a geometric 
representation or as a finite element mesh. 

We will start the implant problem by creating a two-dimensional, deformable 
planar part. We do this by sketching the geometry of the implant. ABAQUS/CAE 
automatically enters the Sketcher when creating the part. ABAQUS/CAE often 
displays a short message in the prompt area indicating what one should do next. 


To create the implant model 


(1) If users did not already start ABAQUS/CAE, type abaqus cae, where abaqus is 
the command used to run ABAQUS. 

(2) Select Create Model Database from the Start Session dialog box that appears. 

(3) In the Module list located under the toolbar (see Fig. 5.38), click Part to enter 
the Part module. 

(4) From the main menu bar, select Part—Create to create a new part. 

(5) Name the part Frame. Choose a two-dimensional, deformable body and a planar 
wire base feature. 

(6) Click Continue to exit the Create Part dialog box. ABAQUS/CAE automatically 
enters the Sketcher. The Sketcher toolbox appears in the left side of the main 
window, and the Sketcher grid appears in the viewport. The Sketcher contains a 
set of basic tools that allow users to sketch the two-dimensional profile of the 
part. 

(7) After we finish the implant geometry in Sketcher, click Done to exit the 
Sketcher. 

The implant part has been created in the screen. 


5.3.2.2 Define Material Properties 


There are three materials to be used in the implant model, which are called Cortical 
bone, Implant, and Trabecul bone, respectively. We use the Property module to 
create materials and to define their properties. Fig. 5.39 shows the Submenus 
available under the Material Manager. 
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To define a material 


(1) In the Module list located under the toolbar, select Property to enter the 
Property module. 

(2) From the main menu bar, select Material > Create to create a new material. The 
Create Material dialog box appears. 

(3) Name the material Trabecul, and click Continue. The material editor appears. (see 
Fig. 5.39) 
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Fig. 5.39 Define Material Properties in ABAQUS/CAE 


(4) Use the menu bar under the browser area of the material editor to reveal menus 
containing all the available material options. When we select a material option, 
the appropriate data entry form appears below the menu. 

(5) Type a value of 1370 for Young's modulus and a value of 0.31 for Poisson's 
Ratio in the respective fields. Use [Tab] or move the cursor to a new cell and 
click to move between cells. 

(6) Click OK to exit the material editor. 

ABAQUS uses section concept to link materials with objects. We need define 
the section properties of a model by creating sections in the Property module. After 
the section creation, we can use one of the following two methods to assign the 
section to the part in the current viewport: 

1) Simply select the region from the part and assign the section to the selected 

region; 

2) Use the Set toolset to create a homogeneous set containing the region and 

assign the section to the set. 

In the implant model, we will create three sections that are behave of Cortical, 
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Implant, and Trabecul, respectively. The sections will be assigned to the three 
relative regions, see Fig. 5.39. 


5.3.2.3 Define the Assembly 


Each part that we create is oriented in its own coordinate system and is independent 
of the other parts in the model. We use the Assembly module to define the 
geometry of the assembly by creating instances of a part and then positioning the 
instances relative to each other in a global coordinate system. Although a model may 
contain many parts, it contains only one assembly. 

We will create a single instance of the implant in this model. ABAQUS/CAE 
positions the instance so that the origin of the sketch that defined the implant 
overlays the origin of the assembly's default coordinate system. 


To define the assembly 


(1) In the Module list located under the toolbar, click Assembly to enter the 

Assembly module. 

(2) From the main menu bar, select Instance Create. The Create Instance dialog box 
ap pears. 
(3) In the dialog box, select Implant and click OK. 

ABAQUS/CAE creates an instance of the implant. In this example the single 
instance of the implant defines the assembly. The implant is displayed in the 1-2 
plane of the global coordinate system (a right-handed, rectangular Cartesian system). 
A triad in the lower-left corner of the viewport indicates the origin and orientation 
of the global coordinate system. 
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Fig. 5.40 Configuring an Analysis in ABAQUS/CAE 
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5.3.2.4 Define an Analysis 


ABAQUS can be used to solve general static and dynamic response problems. In 
this example we are interested in the static response of the implant to both 
horizontal and vertical loads applied at the top center of the implant, with the 
bottom end fully constrained (see Fig. 5.40). This is a stngle event, so only a single 
analysis step 1s needed for the simulation. Thus the analysis will consist of two 
steps overall: 

1) An initial step, in which we will apply boundary conditions that constrain the 

bottom end of the implant; 

2) An analysis step, in which we will apply two concentrated loads at the top 

center of the implant. 

ABAQUS/CAE generates the initial step automatically, but we must use the 
Step module to create the analysis step ourselves. The Step module also allows us 
to request output for any step in the analysis. 

There are two kinds of analysis steps in ABAQUS: general analysis steps, 
which can be used to analyze linear or nonlinear response, and linear perturbation 
steps, which can be used only to analyze linear problems. In this simulation we will 
define a Static, general step. 


To create a static analysis step 


(1) In the Module list located under the toolbar, click Step to enter the Step module. 

(2) From the main menu bar, select Step—>Create to create a step. The Create Step 
dialog box appears with a list of all the general procedures and a default step 
name of Step-1. 

(3) Change the step name to Load Case_1. 

(4) Select “Static, general” as the Procedure type, then continue. The Edit Step 
dialog box appears with the default settings for a static step. 

(5) Click the tabs, Basic/Incrementation/Other, to see its contents; The default 
values provided for the step are acceptable. 

(6) Click OK to create the step and to exit the Edit Step dialog box. 


To request data output 


Finite element analyses can create large amounts of output. ABAQUS allows users 
to control and manage this output so that only data required to interpret the results 
of the simulation are produced. Four types of output are available: 

1) Results stored in a neutral binary file used by ABAQUS/Viewer for 
postprocessing(This file is called the ABAQUS output database file and has 
the extension .odb); 

2) Printed tables of results, written to the ABAQUS data (.dat) file; 

3) Restart data used to continue the analysis, written to the ABAQUS restart 
(res) file; 

4) Results stored in binary files for subsequent postprocessing with third-party 
software, written to the ABAQUS results (.fil) file. 
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By default, ABAQUS/CAE writes the results of the analysis to the output 
database (.odb) file. When a step is created, ABAQUS/CAE generates a default 
output request for the step. We do not need to do anything to accept these defaults. 
The Field Output Requests Manager can be used to request output of variables that 
should be written at relatively low frequencies to the output database from the 
entire model or from a large portion of the model. The History Output Requests 
Manager can be used to request output of variables that should be written to the 
output database at a high frequency from a small portion of the model, for example, 
the displacement of a single node. 

Results can also be written in tabular form to the ABAQUS data (.dat) file. By 
default, no printout is written to the ABAQUS data file. In this implant model, we 
are interested in the displacements (output variable U), the reaction forces at the 
constrainted locations (output variable RF), and the stress in the members (output 
variable S). 


5.3.2.5 Apply Boundary Conditions and Loads to the Model 


In ABAQUS analysis, prescribed conditions such as loads and boundary conditions, 
are step dependent, which means that users must specify the step or steps in which 
they become active. Now that we have defined the steps in the analysis, we can use 
the Load module to define prescribed conditions. 

In structural analysis, boundary conditions are applied to those regions of the 
model where the displacements and/or rotations are known. Such regions may be 
constrained to remain fixed (have zero displacement and/or rotation) during the 
simulation or may have specified, nonzero displacements and/or rotations. The 
directions in which motion is possible are called Degrees of Freedom (DOF). 

In this model the bottom edge portion of the implant is constrainted completely, 
and thus, cannot move in any direction. Fig. 5.41 shows the boundary condition. 
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Fig. 5.41 Define Boundary Condition and Load Case in ABAQUS/CAE 
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To apply boundary conditions to the implant 


(1) In the Module list located under the toolbar, click Load to enter the Load 
module. 

(2) From the main menu bar, select BC-—> Create. The Create Boundary Condition 
dialog box appears. 

(3) In the Create Boundary Condition dialog box: 

1) Name the boundary condition Fixed. 

2) From the list of steps, select Initial as the step in which the boundary 
condition will be activated. All the mechanical boundary conditions specified 
in the Initial step must have zero magnitudes. This condition is enforced 
automatically by ABAQUS/CAE. 

3) In the Category list, accept Mechanical as the default category selection. 

4) In the Type for selected step list, select Displacement/Rotation, and click 
Continue. 

(4) In the viewport, select the vertex at the bottom edge of the implant as the region 
to which the boundary condition will be applied. 

(5) Click mouse button 2 in the viewport or click Done in the prompt area to 
indicate that we have finished selecting regions. 

Then the Edit Boundary Condition dialog box appears. When we are defining a 
boundary condition in the initial step, all available DOFs are unconstrainted by 
default. In the dialog box: 

1) Toggle on U1 and U2 since all translational DOFs need to be constrainted. 

2) Click OK to create the boundary condition and to close the dialog box. 
ABAQUS/CAE displays two arrowheads at the vertex to indicate the 
constrainted DOFs. 


To apply a load to the implant 


After we have constrainted the implant, we can apply two concentrated loads to the 
top center of the implant. In ABAQUS the term “load” (as in the Load module in 
ABAQUS/CAE) generally refers to anything that induces a change in the response 
of a structure from its initial state, including: 

1) concentrated forces; 

2) pressures; 

3) nonzero boundary conditions; 

4) body loads; 

5) temperature (with thermal expansion of the material defined). 

In this simulation two concentrated forces of 100 N are applied in the negative 
l-direction and the negative 2-direction to the top center of the implant, 
respectively (see Fig. 5.41). These loads are applied during the static step we 
created in the Step module. In reality there is no such thing as a concentrated, or 
point, load; the load will always be applied over some finite area. However, if the 
area being loaded is small enough, it is an appropriate idealisation to treat the load 
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as a concentrated load. 
To apply a concentrated force to the implant 


(1) From the main menu bar, select Load-> Manager. The Load Manager appears. 

(2) At the bottom of the Load Manager, click Create. The Create Load dialog box 
ap pears. 

(3) In the Create Load dialog box: 

1) Name the load Force; 

2) From the list of steps, select Apply load as the step in which the load will be 
ap p lied; 

3) In the Category list, accept Mechanical as the default category selection; 

4) In the Type for selected step list, accept the default selection of Concentrated 
force; 

5) Click Continue. 

(4) In the viewport, select the vertex at the top center of the implant as the region 
where the load will be applied. 

(5) Click mouse button 2 in the viewport or click Done in the prompt area to 
indicate that we have finished selecting regions. The Edit Load dialog box 
appears. In the dialog box: 

1) Enter a magnitude of -100.0 for CF]; 

2) Click OK to create the load and to close the dialog box. ABAQUS/CAE 
displays a downward-pointing arrow at the vertex to indicate that the load is 
applied in the negative 1-direction. 

(6) Repeat the above procedure to apply a concentrated load of 100.0 N in the 
negative 2-direction. 

5.3.2.6 Mesh the Model 


We use the Mesh module to generate the finite element mesh. We can choose the 
meshing technique that ABAQUS/CAE will use to create the mesh, the element 
shape, and the element type. The default meshing technique assigned to the model is 
indicated by the color of the model that is displayed when one enters the Mesh 
module; if ABAQUS/CAE displays the model in orange, it cannot be meshed 
without additional assistance. 

In this section we will assign a particular ABAQUS element type to the model. 
Two-dimensional “Quadrilateral” elements will be used to model the implant. These 
elements are chosen because they are ideal for modelling 2D plane strain problem 
such as this implant. 


To assign an ABAQUS element type 


(1) In the Module list located under the toolbar, click Mesh to enter the Mesh 
module. 

(2) From the main menu bar, select Mesh->Element Type. 

(3) In the viewport, select the entire implant as the region to be assigned an element 
type. In the prompt area, click Done when we finish. The Element Type dialog 
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box appears. 

(4) In the dialog box, select the following: 

1) Standard as the Element Library selection (the default). 
2) Linear as the Geometric Order (the default). 
3) Plane strain as the Family of elements. 

(5) In the Line tabbed page, accept the default element type. A description of the 
element type CPE4 appears at the bottom of the dialog box. ABAQUS/CAE will 
now associate CPE4 elements with the elements in the mesh. 

(6) Click OK to assign the element type and to close the dialog box. 

(7) In the prompt area, click Done to end the procedure. 


Create the mesh 


Basic meshing 1s a two-stage operation: first we seed the edges of the part instance, 
and then we mesh it. We select the number of seeds based on the desired element 
size or on the number of elements that we want along an edge, and ABAQUS/CAE 
places the nodes of the mesh at the seeds whenever possible. In this implant model, 
a high mesh density is required in the particular interfaces as shown in Fig. 5. 42. 
Along the gear edge of the implant in Fig, 5. 42, the high mesh density can be 
generated using “seeds along edge” technology in ABAQUS/CAE. 
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Fig. 5.42 Global Mesh and Detailed Mesh of the Implant Model 


To seed and mesh the model 


(1) From the main menu bar, select Seed->Instance to seed the part instance. The 

prompt area displays the default element size that ABAQUS/CAE will use to 
seed the part instance. This default element size is based on the size of the part 
instance. A relatively large seed value will be used so that only corase element 
density will be created in global region. 
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(2) In the prompt area, specify an element size , and press Enter or click mouse 
button 2 in the viewport. 

(3) Click mouse button 2 in the viewport or Done in the prompt area to accept the 
seeding. 

(4) From the main menu bar, select Mesh Instance to mesh the part instance. 

(5) From the buttons in the prompt area, click Yes to confirm that we want to mesh 
the part instance. 


5.3.2.7 Create an Analysis Job and Use Keywords Editor 


We have completely configured the analysis in the above sections from 2.1 to 2. 6. 
In this section, we will move to the Job module to create a job that is associated 
with our model and to add printed output requests using the Keywords Editor. 


To create an analysis job 


(1) In the Module list located under the toolbar, click Job to enter the Job module. 

(2) From the main menu bar, select Job—> Manager. The Job Manager appears. When 
we finish defining our job, the Job Manager will display a list of all jobs, the 
model associated with each job, the type of analysis, and the status of the job. 

(3) In the Job Manager, click Create. The Create Job dialog box appears with a list 
of the models in the model database, see Fig. 5.43. 
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Fig. 5.43 Job Manager to Create a Job for the Implant Model 


(4) Name the job Examg201, and click Continue. The Edit Job dialog box appears. 

(5) In the Description field, type Load Case_1. 

(6) In the Submission tabbed page, select “Data check” as the Job Type to execute 
a check which find out any possible error in the model definition. Click the 
General tab, and toggle all Preprocessor Printout options. Click OK to accept all 
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other default job settings in the job editor and to close the dialog box. 
To add printed output requests using the keywords editor 


(1) From the main menu bar, select Model-> Edit Keywords and select the current 
model. The Edit Keywords dialog box appears containing the input file that has 
been generated for our model. 

(2) Only text blocks with a white background can be edited. Use the scroll bar on 
the right side of the dialog box to find the text blocks where the * EL PRINT 
and * NODE PRINT options are located (toward the bottom of the file). In 
each block enter data by first placing the cursor at the end of the line. For each 
option change the frequency at which printed output will be written to | and 
press Enter to create a new line on which we can specify the output variable 
names. The option blocks should resemble the following when we finish: 

* El Print, freq = 1 

S, 

* Node Print, freq = 1 
U, 

RF, 

(3) When we have finished the whole editing, click OK to exit the Keywords Editor. 

5.3.2.8 Run the Analysis 


When the datacheck analysis completes with no error messages, run the analysis 
itself. To do this, edit the job definition and set the Job Type to Continue analysis; 
then, click Submit in the Job Manager to submit our job for analysis. 

One should always perform a datacheck analysis before running a simulation to 
ensure that the model has been defined correctly and to check whether there is 
enough disk space and memory available to complete the analysis. However, it is 
possible to combine the datacheck and analysis phases of the simulation by setting 
the Job Type to “Full analysis”. 

ABAQUSI/CAE allows users to monitor the running procedure during the running 
time. By selecting “Monitor” in Job manager window, the job monitor window 
appears as Fig. 5.44. For multi increments analysis we are allowed to view the half- 
way results before the analysis is finished. By selecting “Results” in Job manager 
window, the results up to the current increment will be shown in the viewport. It is 
convenient to check if the analysis is going to the right direction before it takes a 
substantial amount of time to go to a wrong direction. 

If a simulation is expected to take a substantial amount of time, it may be 
convenient to run it in a batch queue by selecting Queue as the Run Mode. 


5.3.2.9 Result Analysis in ABAQUS 
ABAQUS job will generate several files that contain different information about the 
job analysis and the job results. 


Three types of message file are available: 
1) Log file is written to record the job procedure; 
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Fig. 5.44 Job Monitor to Monitor Running Analysis 


2) Status file is used to monitor the analysis; 
3) Message file contains detailed information about the simulation’s progress. 


Four types of output are available: 

1) Neutral binary output can be written to the output database (.odb) file, using 
the * OUTPUT option and related suboptions; 

2) Printed output can be written to the printed output (.dat) file. This is only 
available for ABAQUS/Standard; 

3) Restart output can be written to the restart (res) file using the * , RESTART 
option for the purpose of conducting restart analyses; 

4) Results (fil) file output can be written for third-party postprocessing. 


5.3.3 Job Information Files 


5.3.3.1 Log File 


The job log file, examg201.log, contains the messages similar to those shown 
below: 
ABAQUS JOB examg201 
Begin Solver Input File Processor 
06/21/01 15:43:06 
Run C:\ ABAQUS\ 6.2-1\ exec\ pre.exe 
ABAQUS/STANDARD is running on a Category B machine 
and has checked out the Enabling Token 
from the license server on pc001 
06/21/01 15:43:15 
End Solver Input File Processor 
Begin ABAQUSStandard Analysis 
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06/21/01 15:43:15 
Run C:\ ABAQUS\ 6.2-1\ exec\ standard.exe 


ABAQUS/STANDARD is running on a Category B machine 


and has checked out the Enabling Token 
from the license server on pc001 
06/21/01 15:43:41 

End ABAQUS/Standard Analysis 
ABAQUS JOB examg201 COMPLETED 


5.3.3.2 Status File 


Monitor the simulation while it is running by looking at the status file, examg201. 
sta. When ABAQUS has finished the analysis, the status file will contain 


information similar to the following: 


SUMMARY OF JOB INFORMATION: 


IF 


STEP INC ATT SEVERE EQUIL TOTAL TOTAL STEP INC OF DOF 
DISCON ITERS ITERS ‘TIME’ TIME/LPF TIME/LPF MONITOR RIKS 
ITERS FREQ 

l » 4 0 1 1 1.00 1.00 1.000 


THE ANALYSIS HAS COMPLETED SUCCESSFULLY 
5.3.3.3 Message File 


The message file, examg201.msg, contains detailed information about the analysis. In 
the first part of the message file, we can look at the information for the increment 


control parameters. 


STEP 1 INCREMENT l 


STEP 1 


STEP TIME 0.00 
STATICANALYSIS 


This load case is the default load case that always appears 


AUTOMATIC TIME CONTROL WITH— 
A SUGGESTED INITIAL TIME INCREMENT OF 
AND A TOTAL TIME PERIOD OF 
THE MINIMUM TIME INCREMENT ALLOWED IS 
THE MAXIMUM TIME INCREMENT ALLOWED IS 
CONVERGENCE TOLERANCE PARAMETERS FOR FORCE 
CRITERION FOR RESIDUAL FORCE FOR A NONLINEAR PROBLEM 
CRITERION FOR DISP. CORRECTION IN A NONLINEAR PROBLEM 
INITIAL VALUE OF TIME AVERAGE FORCE 
AVERAGE FORCE IS TIME AVERAGE FORCE 


1.00 
1.00 
1.000E-05 
1.00 


5.000E-03 
1.000E-02 
1.000E-02 


ALTERNATE CRIT. FOR RESIDUAL FORCE FOR A NONLINEAR PROBLEM 


CRITERION FOR ZERO FORCE RELATIVE TO TIME AVRG. FORCE 
CRITERION FOR RESIDUAL FORCE WHEN THERE IS ZERO FLUX 
CRITERION FOR DISP. CORRECTION WHEN THERE IS ZERO FLUX 
CRITERION FOR RESIDUAL FORCE FOR A LINEAR INCREMENT 
FIELD CONVERSION RATIO 


2.000E-02 
1.000E-05 
1.000E-05 
1.000E-03 
1.000E-08 
1.00 


VOLUMETRIC STRAIN COMPATIBILITY TOLERANCE FOR HYBRID SOLIDS 


AXIAL STRAIN COMPATIBILITY TOLERANCE FOR HYBRID BEAMS 


1.000E-05 
1.000E-05 
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TRANS. SHEAR STRAIN COMPATIBILITY TOLERANCE FOR HYBRID BEAMS 
1.000E-05 
SOFT CONTACT CONSTRAINT COMPATIBILITY TOLERANCE FOR P> PO 
5.000E-03 
SOFT CONTACT CONSTRAINT COMPATIBILITY TOLERANCE FOR P = 0.0 
0.100 
DISPLACEMENT COMPATIBILITY TOLERANCE FOR DCOUP ELEMENTS 
1.000E-05 
ROTATION COMPATIBILITY TOLERANCE FOR DCOUP ELEMENTS 1.000E-05 


TIME INCREMENTATION CONTROL PARAMETERS 


FIRST EQUILIBRIUM ITERATION FOR CONSECUTIVE DIVERGENCE CHECK 4 
EQUILIBRIUM ITERATION AT WHICH LOG. CONVERGENCE RATE CHECK BEGINS 
8 


EQUILIBRIUM ITERATION AFTER WHICH ALTERNATE RESIDUAL IS USED 9 
MAXIMUM EQUILIBRIUM ITERATIONS ALLOWED 16 
EQUILIBRIUM ITERATION COUNT FOR CUT-BACK IN NEXT INCREMENT 10 
MAXIMUM EQUILIB. ITERS IN TWO INCREMENTS FOR TIME INCREMENT 
INCREASE 4 
MAXIMUM ITERATIONS FOR SEVERE DISCONTINUITIES 12 
MAXIMUM CUT-BACKS ALLOWED IN AN INCREMENT 5 
MAXIMUM DISCON. ITERS IN TWO INCREMENTS FOR TIME INCREMENT 
INCREASE 6 
CUT-BACK FACTOR AFTER DIVERGENCE 0.2500 
CUT-BACK FACTOR FOR TOO SLOW CONVERGENCE 0.5000 
CUT-BACK FACTOR AFTER TOO MANY EQUILIBRIUM ITERATIONS 
0.7500 
CUT-BACK FACTOR AFTER TOO MANY SEVERE DISCONTINUITY ITERATIONS 
0.2500 
CUT-BACK FACTOR AFTER PROBLEMS IN ELEMENT ASSEMBLY 
0.2500 
INCREASE FACTOR AFTER TWO INCREMENTS THAT CONVERGE QUICKLY 
1.500 
MAX. TIME INCREMENT INCREASE FACTOR ALLOWED 1.500 
MAX. TIME INCREMENT INCREASE FACTOR ALLOWED (DYNAMICS) 
1.250 
MAX. TIME INCREMENT INCREASE FACTOR ALLOWED (DIFFUSION) 
2.000 
MINIMUM TIME INCREMENT RATIO FOR EXTRAPOLATION TO OCCUR 
0.1000 
MAX. RATIO OF TIME INCREMENT TO STABILITY LIMIT 1.000 
FRACTION OF STABILITY LIMIT FOR NEW TIME INCREMENT 0.9500 
PRINT OF INCREMENT NUMBER, TIME, ETC., EVERY 1 INCREMENTS 
THE MAXIMUM NUMBER OF INCREMENTS IN THIS STEP IS 100 
LINEAR EXTRAPOLATION WILL BE USED 
CHARACTERISTIC ELEMENT LENGTH 0.297 


PRINT OF INCREMENT NUMBER, TIME, ETC., TO THE MESSAGE FILE EVERY 1 
INCREMENTS 
EQUATION ARE BEING REORDERED TO MINIMIZE WAVEFRONT 
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In the second part of the message file, we can look at the model's behaviour as the following: 


INCREMENT 1 STARTS. ATTEMPT NUMBER 1, TIME INCREMENT 1.00 
EQUILIBRIUM ITERATION it 


AVERAGE FORCE 3.46 TIME AVG. FORCE 3.46 
LARGEST RESIDUAL FORCE 8.664E-11 AT NODE 2399 = DOF | 
LARGEST INCREMENT OF DISP. - 0.194 AT NODE 2240 DOF 1 
LARGEST CORRECTION TO DISP. ~ 0.194 AT NODE 2240 DOF 1 


THE FORCE EQUILIBRIUM RESPONSE WAS LINEAR IN THIS INCREMENT 


ITERATION SUMMARY FOR THE INCREMENT: 

1 TOTAL ITERATIONS, OF WHICH 

0 ARE SEVERE DISCONTINUITY ITERATIONS AND 1 ARE EQUILIBRIUM 
ITERATIONS. 

TIME INCREMENT COMPLETED 1.00, FRACTION OF STEP COMPLETED 1.00 

STEP TIME COMPLETED — 1.00, TOTAL TIME COMPLETED __ 1.00 
THE ANALYSIS HAS BEEN COMPLETED 

At the end of this file the analysis summary is given as below: 


ANALYSIS SUMMARY 


TOTAL OF 1 INCREMENTS 
0 CUTBACKS IN AUTOMATIC INCREMENTATION 
| ITERATIONS 
| PASSES THROUGH THE EQUATION SOLVER OF WHICH 
1 INVOLVE MATRIX DECOMPOSITION, INCLUDING 
0 DECOM POSITION(S) OF THE MASS MATRIX 
1 REORDERING OF EQUATIONS TO MINIMIZE WAVEFRONT 
0 ADDITIONAL RESIDUAL EVALUATIONS FOR LINE SEARCHES 
0 ADDITIONAL OPERATOR EVALUATIONS FOR LINE SEARCHES 
25 WARNING MESSAGES DURING USER INPUT PROCESSING 
0 WARNING MESSAGES DURING ANALYSIS 
0 ANALYSIS WARNINGS ARE NUMERICAL PROBLEM MESSAGES 
0 ANALYSIS WARNINGS ARE NEGATIVE EIGENVALUE MESSAGES 
0 ERROR MESSAGES 
THE SPARSE SOLVER HAS BEEN USED FOR THIS ANALYSIS. 

JOB TIME SUMMARY 

USER TIME (SEC)= 5.8000 

SYSTEM TIME (SEC)= 0.8000 

TOTAL CPU TIME (SEC)= 6.6000 

WALLCLOCK TIME (SEC)= 24 


5.3.4 Job Result Files 


5.3.4.1 Postprocess with ABAQUS/CAE 


Graphical postprocessing is important because of the great volume of data created 
during a simulation. In any realistic model it is impractical for users to interpret 
results in the tabular form of the data file. ABAQUS/CAE allows users to view the 
results graphically using a variety of methods, including deformed shape plots, 
contour plots, vector plots, animations, and X-Y plots. 
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When the job completes successfully, we are ready to view the results of the 
analysis with ABAQUS/CAE. From the buttons on the right edge of the Job 
Manager, click Results. ABAQUS/CAE loads the Visualization module, opens the 
output database created by the job, and displays a fast plot of the model. The fast 
plot is a basic representation of the undeformed model shape and is an indication 
that we have opened the desired file. Alternatively, we can click Visualization in the 
Module list located under the toolbar, select File~>Open, select examg201.odb from 
the list of available output database files, and click OK. 


Viewport: 1 ODE: C: implantlexamg201.0db Viewport: 3) ODE: C:Amplantieaxcang?01.0odb 


S35335335388 


£ 


Fig. 5.45 View Analysis Results 


Fig. 5.45 displays the von Mises stress contour on the deformed model shape. 
From the main menu bar, select Contour plot, using default setting option. 


To compare the undeformed model shape with the deformed model shape 


(1) In the Basic tabbed page of the Deformed Shape Plot Options dialog box, toggle 
on Superimpose undeformed plot; 

(2) Click OK. By default, ABAQUS/CAE plots the undeformed model shape in 
green and the deformed model shape in white. The plot is shown in the right 
part of Fig. 5.45. 


5.3.4.2 Print File Data File 


After the analysis is completed, the data file, examg201.dat, will contain the tables 
of results requested with the * NODE PRINT and * EL PRINT options. The 
results from the implant simulation are shown below. The element output consists 
of a table of axial stress in the elements, while the node output consists of tables of 
nodal displacements and reaction forces. 

In the first part of the print file, ABAQUS writes input deck processing 
information. We can look at the job scale in the print file as below: 
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PROBLEMSIZE 
NUMBER OF ELEMENTS IS 2854 
NUMBER OF NODES IS 2945 
NUMBER OF NODES DEFINED BY USERS 2945 
TOTAL NUMBER OF VARIABLES IN THE MODEL 5890 


(DEGREES OF FREEDOM PLUS ANY LAGRANGE MULTIPLIER VARIABLES) 


If we set output requests using * El print to print the element stress at each 
integration point, we will have element stress output in examg201.dat as below: 


ELEMENTOUTPUT 


THE FOLLOWING TABLE [IS PRINTED FOR ALL ELEMENTS WITH TYPE CPE4 AT THE 
INTEGRATION POINTS 


PT FOOT- 

ELEM ENT NOTE SI] $22 $33 S12 
949 1 3.924 2.216 1.904 — 1.495 
949 2 3.834 2.306 1.904 = 917 
949 3 3.575 2.565 1.904 = A520 
949 4 3.413 24a) 1.904 - 1.841 
950 l 379 2.308 1.887 — 2.914 
950 Z 3.615 2.472 1.887 — 2.883 
950 3 3.1773 21g 1.887 = 2.613 
950 4 3.532 2.000 1.887 - 2.612 
951 291135 4.048 2.096 ~ 3.507 
951 Z 2.969 8:793 2.096 3.332 
1798 3 ~ 6.4473E — 04 1.6320E — 03 0.000 TIDE = 04 
1798 4 — 6.3243E — 04 1.6197E — 03 0.000 9.9955E — 04 
1799 ] — 6.2963E — 04 1.S5747E — 03 0.000 9.3184E — 04 
1799 2 ~ 6.2656E — 04 LS/iGE = 03 0.000 9.3865E — 04 
1799 3 — 6.2530E — 04 1.5703E - 03 0.000 9.2794E — 04 
1799 4 — 6.2228E — 04 1.5673E — 03 0.000 9.3454E — 04 

MAXIMUM 2.8429E — 03 5.4068E — 03 0.000 1.4924E — 03 

ELEMENT 1558 1788 1426 1789 

MINIMUM = 5.15308 — 03 — 7.0632E — 03 0.000 — 6.7281E — 03 

ELEMENT 1437 1558 1426 1558 


We have set output requests using * Node print to print the node 
displacements in examg201.dat. Below is part of the node displacements output in 
the file. 


NODEOUTPUT 
THE FOLLOWING TABLE IS PRINTED FOR ALL NODES 
NODE FOOT-NOTE Ul WZ 
I 2.6145E — 04 — 8.1261E — 04 
2 7.6767E — 05 — 1.4798E — 04 


2 2.4300E — 04 — 5.9756E — 04 
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4 2.0798E — 04 ~ 4,7476E — 04 

5 1.7763E — 04 — 3.8173E ~ 04 

6 L.S107E — 04 — 3.1434E — 04 

7 1.2465E — 04 ~ 2.5713E — 04 

8 9.9932E — 05 ~ 2,.0210E — 04 

9 3.6770E — 05 — 4.0049E — 05 

10 6.5746E — 05 ~ 1,1142E - 04 

ll 5.8316E — 05 — 8.4284E — 05 

12 4 8597E — 05 ~ 6.1147E - 05 

13 2.6480E — 05 ~ 8.8901E — 06 

4569 ~ 4,2065E — 03 — 1.3289E — 03 

4570 — 1.1506E — 02 1.7034E — 03 

4571 — 4.8437E — 03 ~ 1.7563E - 03 

4572 — 1.7639E — 02 — 9,2328E ~ 03 

MAXIMUM 5.7242E — 04 3.2395E — 02 
AT NODE 4557 2120 

MINIMUM — 0.1938 — 3.5224E ~ 02 
AT NODE 2240 1869 


5.3.5 Conclusion 


As an advanced finite element software for the general engineering purpose, 
ABAQUS package has been used widely in all engineering areas to solve problems 
from relatively simple linear analysis to the most challenging dynamic nonlinear 
simulations. To apply ABAQUS on the implant analysis and optimal design, we can 
look at the pontienal areas as the following: 

(1) Linear/nonlinear static analysis; 

(2) Linear/nonlinear dynamic analysis and impact analysis; 

(3) Nonlinear contact analysis between the implant interfaces. 

ABAQUS/CAE can be used to create a complete ABAQUS analysis model. The 
solver (ABAQUS/Standard or ABAQUS/Explicit) reads the input file generated by 
ABAQUS/CAE, performs the analysis, sends information to ABAQUS/CAE to 
allow users to monitor the progress of the job, and generates an output database. 
ABAQUS/CAE is used to read the output database and view the results of the 
analysis. | 

We can perform a datacheck analysis once we have created a model. Error and 
warning messages are printed to the data file. After a successful datacheck analysis, 
estimates of the computer resources required for the simulation are printed to the 
data file. 

Use ABAQUS/CAE to verify the model geometry and boundary conditions 
graphically, using the output database file generated during the datacheck phase. 

It is often easiest to check for mistakes in material properties in the data (.dat) 
file; geometry, loads, and boundary conditions are more easily checked with a 
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graphical postprocessor such as ABAQUS/CAE. 

Tables of results are printed in the data file. The output database file contains 
results for graphical postprocessing with ABAQUS/CAE. 

Always check that the results satisfy basic engineering principles, such as 
equilibrium. 

ABAQUS/CAE allows users to visualize analysis results graphically in a variety 
of ways. 
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